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Engineering  design  requires  increasing  accuracy  in  the 
quantitative  description  of  the  thermodynamic  properties  of 
fluid  mixtures.    To  meet  this  demand,  present  methods  of 
property  calculations  are  more  theoretical  than  empirical; 
they  relate  microscopic  behavior  and  intermolecular  poten- 
tials to  macroscopic  properties  using  statistical  mechanical 
formulations.     Fluctuation  solution  theory  expresses  the 
solution  compressibility  and  component  molar  volumes  and 
activity  coefficients  in  terms  of  integrals  of  the  direct 
correlation  function.     These  are  independent  of  pairwise 
additivity  and  generally  insensitive  to  the  details  of  the 
intermolecular  forces . 

The  present  work  applies  this  theory  to  gas  solubility 
in  mixed  solvents  and  to  liquid-liquid  equilibria.     It  in- 
cludes an  extensive  collection  of  experimental  direct 


xvii 


correlation  function  integrals  for  liquids  and  makes  some 
suggestions  for  improving  the  original  model.     In  particu- 
lar,  ideal  solution  direct  correlation  function  integrals 
have  a  simple  relation  to  composition.     In  nonideal  mixtures, 
while  the  form  may  appear  more  complex,   simple  relations 
are  also  consistent  with  the  data. 

It  has  been  found  that  the  fluctuation  solution  form 
of  the  Van  der  Waals  model  describes  liquid-liquid  equi- 
libria in  binary  systems  as  well  as  any  other  liquid  model 
currently  in  use.     It  is  also  excellent  for  gases  in  mixed 
solvents.     However,  its  use  of  only  one  parameter  is  a 
severe  handicap  when  the  model  is  applied  to  multicomponent 
liquid-liquid  systems.     To  improve  the  situation,  different 
alternatives  have  been  examined  but  none  has  proven  entirely 
satisfactory.     Various  choices  of  hard  sphere  expression 
in  the  perturbed  hard  sphere  model  leads  to  the  same  inade- 
quate results  in  mixtures.     Hard  convex  body  equations  make 
no  improvement.     The  major  problem  is  indicated  by  the 
fact  that  at  infinite  dilution  in  a  binary  system  the 
model  does  not  account  for  the  reality  of  the  unlike  inter- 
actions.    The  introduction  of  a  quadratic  term  in  density 
may  be  a  way  to  solve  this  problem.     A  study  of  gases  in 
mixed  solvents  is  also  included. 
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CHAPTER  1 
INTRODUCTION 


Design  of  chemical  processing  equipments  requires 
accurate  phase  equilibria  descriptions .     To  cover  all 
necessary  ranges  of  conditions  would  demand  a  large  amount 
of  experimental  work.    A  way  to  avoid  the  expensive  and  the 
time-cons\aming  path  of  experiments  is  to  predict  the  behav- 
ior from  a  minimum  of  fundamental  data,  e.g.,  multicompo- 
nent  system  properties  from  pure  component  and  binary 
properties . 

Different  classes  of  fluid  phase  equilibria  arise  from 
the  physical  and  chemical  nature  of  the  mixture.    Thus,  it 
is  possible  to  talk  about  two  or  more  phases  among  gases, 
liquids,  and  solids  including  multifluid  systems.     The  sub- 
stances can  be  categorized  as  electrolytes  and  nonelectro- 
lytes  with  important  subdivisions  of  nonpolar,  polar,  and 
aqueous  solutions  of  discrete  or  continuous,  and  small  or 
polymeric  molecules.     Further  characterization  and  division 
are  possible  using  the  conditions  of  temperature  and  pres- 
sure on  the  mixture.     This  work  is  an  attempt  to  build  a 
predictive  method  for  liquid-liquid  equilibria  of  nonelec- 
trolyte  systems  of  nonpolymeric  molecules . 

Thermodynamics  tells  us  that  at  two-phase  equilibrium 
the  chemical  potential  of  each  component  is  the  same  in 
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2 

both  phases.     These  chemical  potentials  depend  on  thermo- 
dynamic variables  of  the  system  such  as  temperature,  pres- 
sure, volume,  and  chemical  composition.     The  problem  is 
to  find  a  good  relation  to  expression  for  the  chemical 
potential  in  terms  of  these  variables.     In  order  to  simplify 
the  problem,  it  is  common  to  introduce  functions  such  as 
the  fugacity  and  activity   (Prausnitz,   1969)  .     Hence,  the 
phase  equilibria  problem  reduces  to  finding  a  model  for 
the  fugacity  or  the  activity. 

The  usual  predictive  methods  in  liquid-liquid  phase 
equilibria  can  be  divided  in  two  classes:     those  that  use 
the  same  equation  for  both  phases  and  those  that  use  dif- 
ferent equations  for  each  phase.     The  first  class  works 
with  fugacity  coefficients  and  writes  them  from  a  P-v-T-x 
equation  of  state.     The  second  class  works  with  activities 
or  activity  coefficients  in  the  liquid  phase  and  writes 
them  from  excess  Gibbs  energy  expressions. 

The  fluctuation  solution  theory  approach  used  here 
does  not  belong  to  either  of  these  two  classes.     It  employs 
for  the  liquid  phase  an  equation  of  state  only  for  the 
liquid  phase  rather  than  both  gas  and  liquid  as  the  first 
class  does,  and  it  is  not  derived  from  an  excess  Gibbs 
energy  model.    Another  distinctive  and  important  point  is 
the  generality  of  the  approach  where  expressions  for  the 
species,  whether  gases,   liquids,  or  ions,  are  evaluated  in 
the  same  way. 
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In  the  fluctuation  solution  theory,  thermodynamic 
properties  are  obtained  by  integration  of  a  model  for 
derivatives  of  the  properties  so  the  approach  is  likely  to 
be  less  sensitive  to  errors  in  modelling.     In  addition,,  the 
theory  employs  temperature  and  component  density  as  in- 
dependent variables,  as  opposed  to  pressure,  temperature  and 
mole  fraction  which  are  normally  the  choice.     This  allows 
the  possibility  of  simpler  expressions  since,  as  will  be 
shown,  the  results  are  more  sensitive  to  the  density. 

The  goal  in  this  work  is  to  build  on  the  success  of 
the  one-parameter  corresponding  states  model  for  pure  com- 
ponent liquids  of  Brelvi   (1973)   and  the  two-parameter  cor- 
responding states  model  for  solutions  containing  super- 
critical components  of  Mathias  (19  78) .     Other  related  work 
has  dealt  with  solutions  of  "reactive  components"  including 
electrolytes   (Perry,   1980)   and  groups   (Telotte  and  O'Connell, 
1982)  . 

Chapter  2  describes  the  fluctuation  solution  theory 
and  how  it  can  be  used  to  solve  phase  equilibrium  problems. 
Using  the  grand  canonical  ensemble,  direct  correlation 
function  integrals   (DCFI)   are  related  to  the  concentration 
derivatives  of  the  chemical  potentials  of  the  components . 
When  these  expressions  are  integrated,  equations  for  the 
thermodynamic  properties  of  interest  are  obtained.  The 
chapter  ends  with  the  stability  conditions  for  multicompo- 
nent  systems  in  this  new  theory. 
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Chapter  3  displays  direct  correlation  function  inte- 
grals for  many  ideal  and  nonideal  binary  systems.  These 
integrals  are  obtained  from  volumetric  and  excess  Gibbs 
energy  correlation  of  measured  data,  so  these  integrals 
are  considered  to  be  experimental.     The  apparent  composi- 
tion dependence  in  some  cases  is  more  complex  than  expected. 
However,   it  is  shown  that  a  simpler  form  of  the  DCFI  pro- 
duces essentially  the  same  property  results.    This  indicates 
the  direction  of  sensitivity  for  DCFI  modelling  of  thermo- 
dynamic properties. 

Chapter  4  presents  a  model  for  direct  correlation 
function  integrals.     Mathias   (1978)   used  it  to  describe 
the  solubility  of  gases  in  liquids.     The  model  writes  the 
DCFI  in  the  Van  der  Waals  form,  a  hard  sphere  expression 
plus  a  correction  term  linear  in  the  density.     In  this 
chapter,  we  analyze  two  modifications  to  this  model.  One 
of  them  introduces  a  parameter  in  the  hard  sphere  equation. 
The  correction  term  is  adjusted  to  make  the  model  accurately 
describe  pure  component  properties .     The  other  replaces  the 
hard  sphere  equation  by  a  hard  convex  body  expression. 
Finally,  we  show  application  of  the  Mathias  model  to  the 
solubility  of  gases  in  binary  solvents. 

Chapter  5  describes  application  of  the  model  of 
Chapter  4  to  liquid-liquid  equilibria.     This  is  a  more 
severe  test  to  the  fluctuation  solution  theory  because  the 
temperature  effects  and  the  balance  of  like  and  unlike 
attractive  forces  are  more  important  here  than  in  gas-liquid 
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equilibria.    As  given,  the  model  is  limited  to  one  binary 
parameter.     As  a  result,   the  correlation  of  ternary  systems 
is  poor.     Different  ways  are  tried  to  use  two  binary 
parameters . 

Chapter  6  suggests  a  more  drastic  modification  to  the 
model  than  in  Chapter  5.     First,  we  analyze  the  effective- 
ness of  the  hard  sphere  mixture  expressions  as  a  reference 
for  describing  direct  correlation  function  integrals.  Then, 
we  address  the  problem  of  the  infinite  dilution  behavior 
of  the  direct  correlation  function  integrals  which  in 
binary  systems  depend  on  the  interaction  between  the  com- 
ponents.    Finally,  we  look  at  expressions  for  the  direct 
correlation  function  integrals  that  are  obtained  from  some 
common  excess  Gibbs  energy  models . 


CHAPTER  2 

THERMODYNAMICS  AND  FLUCTUATION  SOLUTION  THEORY 
2 . 1  Introduction 

This  chapter  describes  fluctuation  solution  theory  and 
how  it  may  be  employed  for  solving  phase  equilibria  problems. 

Section  2 . 2  employs  the  grand  canonical  ensemble  in  two 
ways:     first,  to  relate  composition  fluctuations  to  deriva- 
tives of  the  chemical  potential  of  the  components,  and 
second,  to  relate  the  composition  fluctuations  to  integrals 
of  the  pair  distribution  function  of  the  molecular  pairs 
present  in  a  solution.     When  the  composition  fluctuations 
are  eliminated,  direct  correlation  function  integrals  are 
related  to  the  derivatives  of  the  basic  equilibria  equa- 
tions.    Section  2.3  shows  how  these  quantities  are  related 
to  the  thermodynamic  properties  of  interest.     Section  2.4 
shows  how  stability  conditions  are  easily  evaluated  in  the 
new  theory. 

Before  going  into  the  details  of  the  derivation,  we 
note  the  two  most  desirable  features  of  the  theory. 

1.     Thermodynamic  properties  are  obtained  by  integra- 
tion of  a  simple  model  instead  of  the  normal  situation  that 
requires  a  derivative  of  a  complex  expression.  Numerical 
analysis  indicates  that  the  present  approach  is  less  sen- 
sitive to  errors  in  modelling. 

6 


7 

2.     The  fluctuation  solution  theory  employs  tempera- 
ture and  component  density   (component  number  density  at 
microscopic  level)  as  independent  variables.     From  a  statis- 
tical mechanical  point  of  view,  this  set  is  the  most 
appropriate  for  theoretical  work.     In  contrast,  pressure, 
temperature,  and  mole  fraction  are  the  standard  choice  in 
thermodynamics,  but  their  use  causes  difficulty  in  develop- 
ing a  complete  equation  of  state  for  either  the  fluid  phases 
or  the  liquid  alone. 

2.2    Statistical  Thermodynamics 

The  basic  references  for  this  chapter  are  Kirkwood  and 
Buff   (1951),  O'Connell   (1971a),  and  Rowlinson  and  Swinton 
(1982) .     The  interactions  which  determine  the  bulk  proper- 
ties of  the  matter  are  electromagnetic;  they  arise  from  the 
Coulomb  interactions  between  nuclei  and  electrons.     One  way 
to  predict  the  properties  of  liquid  would  be  to  solve  the 
many  body  Schrodinger  equation  describing  the  motion  of  the 
nuclei  and  electrons.     To  facilitate  the  task  there  are 
three  basic  simplifications. 

The  first  is  the  Born-Oppenheimer  approximation  for  a 
static  configuration  of  the  nuclei  relative  to  the  electronic 
motion. 

The  second  approximation  is  to  neglect  any  coupling 
between  vibrational  and  center  of  mass  motions  of  the  mole- 
cule.    For  complex  molecules  this  means  that  they  are  rigid. 
The  position  and  orientation  of  a  rigid  molecule  i  are 
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defined  by  six  coordinates,  the  center  of  mass  coordinates 
r^  and  the  Euler  angles  S^^;  in  the  case  of  axially  symmetric 
molecules,  two  angles  are  sufficient  to  define  the  orienta- 
tion.   As  the  intermolecular  forces  are  dependent  both  on 

position  and  orientation,  we  shall  denote  the  potential 

N  N 

energy  function  of  N  molecules  by  the  symbol  U^{r_  ,  ft  ). 

The  third  approximation  is  that  we  can  use  classical 
mechanics  and  classical  statistical  mechanics  in  describing 
the  behavior  of  the  molecules. 

We  consider  an  open  region  of  volume  V,  which  is  a 
part  of  a  system  of  infinite  extent,  and  we  select  the 
grand  canonical  ensemble  which  is  the  most  convenient  for 
deriving  theoretical  results  to  represent  its  statistical 
behavior.     The  volume  contains  exactly  N^,N2,...,N^ 
molecules  of  the  n-species  of  the  multicomponent  system. 
The  system  is  defined  by  the  variables  temperature  (T) , 
total  volume   (V) ,  and  chemical  potential  of  the  n  components 
.     We  can  write  the  probability  that  the  n-component 

system  will  contain  Nj^,N2  molecules  of  the  n  species 

as 

1 

^  exp   (  j-^     Z  N^y^)   Q(T,V,N^,  .  .  .  ,N^) 
^(N^  ,N-,  .  .  .N  )       "   (2-1) 

where  k  is  the  Boltzmann  constant;  Q  and  H  are  the  canonical 
and  grand  canonical  partition  functions  defined  by 

Q  =        .   .  /  e-«/^^  dp^  dqN  (2-2) 


with  H  as  the  hamiltonian  of  total  energy  of  the  system  or 
the  result  of  adding  the  kinetic  energy  function  T  to  the 
potential  energy  function  U. 


N' 
n 

^    1    E  N.u. 
kT  1    ^  ^ 


Q  (2-3) 


1  n 


From  the  definition  of  the  average  ntamber  of  molecules 
we  can  write 


<N . >  =         E  N.P 

^         N,  .  .  .N  ^ 

1  n 


kT  3E 


(2-4) 


From  the  definition  of  the  average  number  of  pairs  of  mole- 
cules we  can  write 


<N.N .  > 
1  3 


I  N.N.P 


N,  .  .  .N 

1  n 


1  : 


2  2 
(kT)^  9^H 


9<N.  > 

<N  .  >    <N  .  >   +   ^T-ni- 

1         J  3yj/kT 


(2-5) 


Distribution  functions  allow  a  complete  description  of 
the  microscopic  structure  of  liquids  and  give  a  quantitative 
measure  of  the  correlations  between  the  positions  of 
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different  particles.  The  most  important  distribution  func- 
tions are  the  singlet  and  pair  molecular  functions. 

For  a  pure  fluid,  the  probability  density  p  that  M 
particular  molecules  lie  in  the  elements  dr^,  dJ^^;  dr2, 
d^2>  •  •  -  >  ^£m'  ^^hi        given  by 


-     V     exp   (N3y)    r       r  I  £3TT  \ 

p  =     i-     —  J...  J     exp   (-6U  ) 

N=0     N!  H 


There  are  N{N-1) . . . (N-M+1)  different  sets  of  molecules  which 
can  occupy  the  volume  elements  so  that  the  total  probability 
density  that  any  N  molecules  occupy  these  elements  is  given 
by 

(m)  M  (M) 

n'    {r^,g^}...7  r^,q^)  =  P  g    '  {r^,q_^; . . . ;  E^'^^ 

=  Ml   2   (2-7) 

^'      (N-M)  ! 

The  functions  n^^^  and  g^^^  are  called  M-molecular  dis- 
tribution functions.     For  an  isotropic  fluid,  n^"^^   is  just 
the  number  density  and  will  depend  only  on  the  distance 
r^2  ~  \e.i~'£.2^'  direction  and  the  orientations  of  the 

molecules  with  respect  to  each  other.     The  pair  distribu- 
(2) 

tion  g       (£3_/R2^/£2'— 2^   plays  an  important  role  in  the  theory 
of  liquids.     It  goes  to  unity  as  the  distance  between  the 
molecules  becomes  large,  indicating  that  no  correlation 
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exists.     For  spherical  molecules  the  pair  distribution  de- 
pends on  the  distance  only,  and  it  is  called  the  radial 
distribution  function.     For  nonspherical  molecules  we 
define  the  radial  distribution  function  by 


In  equation  (2-8)  the  subscript  s  indicates  an  angle 
average . 

Other  important  functions  are  the  total  correlation 
function,  h  =  g-1,  and  the  direct  correlation  function,  c, 
defined  by  the  Ornstein-Zernike  equation  for  a  pure  homo- 
geneous fluid. 


(The  superscript  2  that  above  had  indicated  a  pair  molecu- 
lar function  has  now  been  dropped.)     We  can  write  the  angle 
averaged  version  as 


gg(r^2)  =  g(r^2'Ri'5l2^^^1^^ 


(2-8) 


+  p  /  h'^^2.'—3'—l'-3^   c  (r2,r2  ,f2_2,fl3)  dr^dfi 


(2-9) 


^s(^12)   =  ^s(^12)   ^  P  ^  hs(^13)^s<^23)'^2:3 


(2-10) 


Figure   (2-1)   compares  the  range  and  shape  of  the  pair 
potential,  the  radial  distribution  function,  and  the  direct 
correlation  function  of  a  spherically  symmetric  molecule. 
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For  a  system  of  n  components,  we  define  functions  con- 
sistent with  those  for  pure  fluids,  the  single  molecular 
distribution  function,  n^^^^(rj^)  as  the  time-averaged 

density  of  molecules  of  type  i  at  the  position  r^^.  The 

(2) 

pair  molecular  distribution  function,  ^^j^j  ^£^'£2^' 
defined  as  the  average  density  of  molecular  pairs  of  type 
i  at  r^  and  type  j  at  Integration  of  the  above  dis- 

tribution functions  over  the  volume  of  the  system  gives, 
with  the  subscript  s  eliminated. 


/  n^^-'-Nrj^)dr^  =  <N^>  (2-11) 


(2) 

//  n. .     '{rT,r,)dr,  dr^  =  <N.N.>  -  6..  <N.>  (2-12) 


The  radial  distribution  function  is  defined  by 


^ij      (£i'£2)  (if,    .     (1)  ,  , 

n^'   ' (r^)nj     ' (£2) 


=  gij(£i'£2)  (2-13) 


Using  this  radial  distribution  function  in  equation  (2-12) , 
we  can  write 

,  <N.N.>  -  <N.><N.>     5. . 

k  ^    (^ij  (^12)  -  1)         =       <N^<N.>  '       ^    -  ^  (2-14) 

ID  1 


Combining  equations   (2-5)   and  (2-14)  we  get 
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,  3<N.> 
1  1 


<N^>  8y  ./kT 


=  p     /   (gi.(r,  J-1)  dr2+6  (2-15) 


As  the  volume  of  the  ensemble  is  constant,  the  above  equa- 
tion can  be  written 


1     '  Pi 


P  3yj/kT 


=  x.x.p/   (g. . (r, -1)  dr- 


+  x^  (2-16) 


It  is  possible  to  work  with  equation  (2-16)   in  develop- 
ing models  for  liquid  solutions,  but  the  approach  leads  to 
complex  equations  which  are  ill-defined  at  the  critical 
point.     Therefore,  it  is  desirable  to  express  equation  (2-16) 
in  terms  of  the  direct  correlation  function  which  has  the 
same  range  as  the  intermolecular  potential  function.  Be- 
sides, there  are  strong  indications  that  the  direct  cor- 
relation functions  are  less  difficult  to  model   (Gubbins  and 
O'Connell,  1974)  .     For  a  system  of  n-components  the  equiva- 
lent equation  to   (2-10)  is 


^ij(^12)   =  5ij(^12)    -  1 

(2-17) 

=  c..(r3L2)   +  Pk^^ik(^13)hjk(^23)^£3 


We  must  emphasize  that  in  equation  (2-17)  we  have  dropped 
the  subscript  s  that  indicates  an  angle  average  and  the 
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superscript  (2)   that  indicates  that  they  are  pair  molecu- 
lar functions.     Integrating  equation  (2-17)  over  all 
gives 


where  C  is  a  matrix  of  direct  correlation  function  integrals 

j .     The  term  H  is  a  matrix  with  total  correlation  func- 
tion integrals,  H^^;  X  is  a  diagonal  matrix  with  elements  X^^ 
equal  to  the  mol  fraction  x^.     Postmultiplying  equation 
(2-19)  by  the  matrix  X  and  adding  I_ 

I+HX=CX     +     CXHX     +     I  (2-20) 

Equation  (2-20)  can  also  be  rearranged  to  give 


p  /  h.  .(r^2)dr2  =  p  /  c.j(r3_2)dr2 


(2-18) 


In  matrix  notation  equation  (2-18)   can  be  written  as 


H  =  C  +  C  X  H 


(2-19) 


(I 


H  X)    (I    -    C  X)      =  I 


(2-21) 


We  can  write  equation   (2-16)   in  matrix  form  as 


D     +     X   (I     +  H  X) 


(2-16a) 
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where  the  elements  of  D  are  the  left-hand  sides  of  equa- 
tion (2-16)  . 

Combining  equations   (2-16a)   and  2-21) ,  we  get 


d"-"-    =    (x"-"-  -  C) 


(2-22) 


After  dividing  by  p,  each  element  of  the  matrix  equation 
(2-22)  is 


N  ^Uj/kT 

P  3 


T,V,N 


k?^j 


ay^/kT 
~3p1 


T,p 


k?^j 


(2-23) 


The  Gibbs-Duhem  equation  at  constant  temperature 
relates  changes  in  pressure  and  chemical  potential  as 


VdP  = 


n 

Z  N.dy. 
1     ^  ^ 


(2-24) 


Dividing  equation  (2-24)  by  dp^  at  fixed  T  and  Pj^^j 


gives 


8P/kT 


Pk?^j 


n 


ay^/kT 

~3p^ 


(2-25) 


T,p 


k?^j 


Finally,  using  equation   (2-23)   in  equation   (2-25)  we  get 
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N  3p/kT 
p     3N . 


T,V,N. 


n 


=  Z  X.    (1  -C.  . 
1     ^  ^3 


(2-26) 


where  x.   is  given  by 


Equations  (2-23)  and  (2-26)   are  the  basic  equations  of 
the  fluctuation  solution  theory  when  T  and  £  are  the  in- 
dependent intensive  variables  of  the  system. 


Prior  to  using  the  basic  equations  it  must  be  rec- 
ognized that 

a)  The  most  common  approach  in  modelling  thermodynamic 
properties  of  solutions  requires  taking  derivatives  of  the 
system  free  energy.     On  the  contrary,  here  we  are  going  to 
use  the  integrals  of  a  model  for  the  direct  correlation 
function  integrals.     In  general  this  is  expected  to  produce 
results  that  are  less  sensitive  to  the  model  approximation. 

b)  The  most  common  representation  of  thermodynamic 
properties  uses  temperature,  pressure,  and  mole  fraction 
of  the  components  of  the  system  as  the  independent  vari- 
ables.    But  equations   (2-23)   and  (2-26)  use  temperature  and 
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concentration  of  the  components.     This  means  that  in 
solving  the  system  of  equations  the  pressure  equation 
(equation  (2-26) )  must  be  used  to  find  the  density  of  the 
solution,  p   (O'Connell,  1981). 

The  isothermal  differential  of  any  thermodynamic 
property  Z,  which  is  a  function  of  temperature  and  com- 
ponent density,  is  given  by 


n 

dZ    =  E 


3Z 


1  ''i 


T,  P    .     y  . 


dp. 


(2-27) 


Thus,  the  isothermal  change  in  the  thermodynamic  property 
Z  from  an  initial  state  r  to  a  final  state  f  is 


Z^  -  Z^ 


f 
Pi 

I  J. 


dp. 

^1 


T,  p  . 


dp. 


(2-28) 


Then  for  the  change  in  pressure,  using  equation  (2-26) ,  we 
obtain 


RT 


(P^-P^) 


Pi  P.C..(T,p) 
-  Z  /         S  -iLJJ  =_ 

i     P^    j  P 


dp. 


(2-29) 


And  for  the  change  in  chemical  potentials,  using  equations 
(2-23) ,  we  get 

f 

dp.  (2-30) 


P. 


Pi  ^ 


T,p 
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In  solving  phase  equilibrium  problems  it  is  more  con- 
venient to  work  with  the  fugacity  or  activity  coefficient 
than  with  chemical  potential.  As 


U°(T,p_^)    +  RT  In  Yj_(T,£) 


(2-31a) 


U.  (T,p   )    +  RT  In   — 

fi(T,p'') 


(2-31b) 


We  can  rewrite  equations  (2-31a)  and  (2-31b) 


3  In 


T,  p 


(2-32) 


9  In  f. 
 1 

dp. 


T,p 


(2-33) 


Thus  for  the  change  in  fugacity  we  get 


In  f  r  -  In  f .     =     In  —  -  Z  r 


dp. 


(2-34) 


And  for  the  change  in  activity  coefficient  we  set 

f 

f             r                             C  (T,p) 
In  Yi  -  In  =     Z  /       [  1  -  ^J—         ]  dp  . 

j      r  P  J 


(2-35) 
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In  summary,  the  procedure  for  finding  pressure  and 
chemical  potentials  (or  their  equivalent  properties , 
fugacities,  or  activity  coefficients)   requires  three 
steps : 

1)  Choose  an  appropriate  lower  integration  limit  or 
reference  state  for  the  set  of  partial  differential  equa- 
tions. 

2)  Solve  iteratively  equation  (2-29)   for  the  final 
total  density. 

3)  Calculate  the  activity  coefficient  for  component  i 
using  equation  (2-35) . 

2.4    Applications  to  Phase  Equilibria  Problems 

The  thermodynamic  equilibrium  between  m  different  phases 
in  an  n  multicomponent  system  is  given  by 

1  2  m  .       ,  /  -     ^ » 

yj_    =  =     .   .   .     =  i  =  1,   .   .   .,  n  (2-36) 

where        is  the  chemical  potential  of  component  i  in  phase 
j.     As  equation  (2-36)   is  not  very  useful  in  the  modeling 
real  systems,  it  is  convenient  to  replace  it  by  an  expres- 
sion that  uses  the  fugacity  as  defined  by 

yj_    =    u°  +  RT  In  -i  (2-37) 

■^i 

where  f^  is  the  fugacity  of  component  i  and  the  superscript 
o  indicates  the  reference  for  measuring  the  thermodynamic 
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properties.    The  reference  states  are  chosen  different 
for  each  phase.     The  only  assiomption  is  that  these  ref- 
erence states  are  at  the  same  temperature.     Then,  using 
equation  (2-37)   in  equation  (2-36) ,  the  equilibria 
between  phases  could  be  described  by 


f^  =  f?  =  .   .   .  =  f^        i  =  1,   .   .   .   ,  n  (2-38) 


Another  possibility  is  working  with  the  activity  as  defined 
as 


a.    =   ^  (2-39) 
i 

or 

Uj^    =   y°  +  RT  In  a^  (2-40) 


We  can  rewrite  equation  (2-38)  in  terms  of  the  activities 
^i 


^ol  ^1  -o2  2 

f.     a.     =     f.  a. 
11  11 


^om  m  .  T  V 
=     f^    a^        1  =  1,   .   .    .    ,  n  (2-41) 


Another  alternative  is  to  introduce  the  activity  coeffi- 
cient defined  by 
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=    YiX.  (2-42) 
With  this  last  definition  equation  (2-38)  can  be  written  as 


.ol     11  ^o2     2  2 

f.      v.    X.      =     f.  Y-X. 
1       '  X      1  X       '  X  X 


=  Yi  X  =  X,   .    .    .   ,  n  (2-43) 


where  the  superscript  o  indicates  the  reference  state,  the 
superscript  j  the  phase,  and  the  subscript  i  the  component. 
Up  to  this  point  we  have  presented  the  basic  equations 
without  any  reference  to  two  basic  problems,  the  reference 
state  and  the  variables  of  which  the  properties  depend  upon. 

The  standard  choice  is  to  adopt  a  reference  state  that 
is  composition  independent.     An  ideal  solution  is  defined 
as  the  basis  for  real  solutions 


v]^    =     y°(T,P)   +  RT  In  X.  (2-44a) 


or 


4^    =     X.   f°  (2-44b) 


where  we  have  chosen  the  variables    temperature,  pressure, 
and  mole  fraction  to  describe  the  thermodynamic  properties. 

It  is  found  experimentally  that  all  dilute  solutions, 
where  the  mole  fractions  of  one  component  are  near  zero, 
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behave  like  ideal  solutions.    These  ideal  mixtures  follow 
Henry's  law.     On  the  other  hand,  there  are  solutions  that 
behave  ideally  near  the  pure  component  composition.  These 
mixtures  follow  Raoult's  law.     Deviations  from  ideal  solu- 
tion are  contained  in  the  activity  coefficient. 

For  the  particular  case  of  the  solubility  of  gases  in 
liquids,  equation  (2-38)   for  component  i  takes  the  following 


where  we  have  chosen  temperature,  pressure,  and  mole 
fraction  as  variables,  y.  being  the  set  of  vapor  mole  frac- 
tions and  X  being  the  set  of  liquid  mole  fractions. 

The  fugacity  of  the  liquid  phase  can  be  expressed  in 
several  alternatives;  the  more  important  ones  are 


form: 


fl{T,P,Y)     =  fi(T,P,x) 


(2-45) 


(i)^(T,P,x)P 


(2-46) 


P 


v^ (T,P,x) 


=    x^  Yj_(T,x)    f°(T)   exp  [/ 


RT 


T,x 


dP] 


(2-47) 


=  Yi   (T,P,x)    f°  (T) 


(2-48) 


Equation  (2-46)   is  used  when  an  equation  of  state  is 
available  to  describe  both  phases.     Equations   (2-47)  and 
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(2-48)  can  be  written  from  equation  (2-43) .     In  equation 
(2-47)  the  exponential  term  is  the  Poynting  correction. 
Equation  (2-48)   is  more  appropriate  for  working  with  the 
equations  written  in  section  2.3,  where  equation  (2-35) 
provides  the  activity  coefficient  necessary  for  equation 
(2-48)  .     The  dependence  of        in  equation  (2-48)  can  be 
taken  as  either  (T,P,x)  or  as  (T,p_)  ,  the  latter  being  the 
natural  variables  for  fluctuation  solution  theory. 

For  the  particular  case  of  the  liquid  immiscibility , 
equation  (2-38)   for  component  i  takes  the  following  form: 

f^   (T,P,x')     =     f^   (T,P,x")  (2-49) 

where  we  have  again  chosen  temperature,  pressure,  and  mole 
fraction  as  variables,  x  is  the  set  of  liquid  mole  frac- 
tions, and  the  superscripts  '  and  "  indicate  two  different 
liquid  phases.     If  the  same  reference  state  is  chosen  for 
both  phases,  equation  (2-49)  can  be  written  as 

(T,P,x')   x^     =     Yi    (T,P,x")    x^  (2-50) 

Equation  (2-50)  can  be  used  with  equation  (2-35)  of  section 
2.3. 

Comparison  of  equations   (2-50)   and  (2-45)   shows  why 
the  problem  of  liquid-liquid  equilibria  is  more  difficult. 
First,  there  is  no  independent  determination  of  the 
component  fugacity  as  in  the  vapor-liquid  case  with  the 
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vapor  fugacity.     The  same  fact  does  not  allow  us  to  check 
the  consistency  of  data.     Second,  pure  component  fugacity 
values  have  no  influence,  and  all  variations  arise  from 
nonideality,  expressed  by  the  activity  coefficient.  This 
increases  sensitivity  of  the  results  to  the  model. 

Equation  (2-36)  and  the  others  derived  from  it  say 
nothing  about  equilibrium  stability.     Any  thermodynamic 
density  variable  (entropy,  internal  energy,  Helmholtz  free 
energy,  etc.)  may  be  employed  to  express  stability  criteria. 
The  Helmholtz  free  energy  is  chosen  here  primarily  because 
its  natural  variables  agree  with  variables  in  which  direct 
correlation  function  integrals  are  written.     We  can  state 
that  a  system  is  in  a  stable  equilibrium  state  at  constant 
volume,  temperature,  and  number  of  moles  when  its  Helmholtz 
free  energy  is  at  a  minimum.     Or,  in  other  words,  at 
equilibrium  any  small  perturbation  at  fixed  T,V,N  results 
in  an  increase  in  the  total  Helmholtz  free  energy.  For 
such  small  perturbation,  A  can  be  expanded  in  a  Taylor 
series 

AA    =     5A  +  6^A  +  S^A  +  .    .    .   >  0  (2-51a) 
For  equilibrium  to  exist, 
5A    =  0 

and  for  the  system  to  be  stable. 
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6^A  > 


where  j  is  the  first  nonvanishing  variation.     Choosing  j  =2, 


6^A  >  0  (2-52) 


We  mentally  divide  the  system  in  two  subsystems,  one 
of  them  very  small.     We  generate  a  perturbation  by  allowing 
a  virtual  transfer  of  volume  and  masses  across  the  hypo- 
thetical surface  separating  the  two  subsystems.     As  one  of 
the  subsystems  is  very  small  and  both  are  at  equilibrium, 
equation  (2-51a)   transforms  to 


AA  =  6A 


=  N  da  >  0 


(2-51b) 


where  a  is  the  molar  free  Helmoltz  energy  (small  letters 

represent  molar  properties) . 

For  convenience  we  prefer  to  work  with  extensive 

properties.     Introducing  the  total  number  of  moles  inside 
2 

d  a,  we  get  the  second  differential  of  the  free  Helmoltz 
energy  which  depends,  here,  on  volume  and  (n-1)  niimber  of 
moles   (Modell  and  Reid,  1983) ,  where 

2  T 

d  A    =    P     D  P  (2-53) 

where 
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dV 


P  = 


dN 


n-1 


(2-54) 


D  = 


2 


9V3N 


n-1 


8N     T  9V 
n-1 


^^n-l^Vl 


(2-55) 


The  stability  criteria  are  the  criteria  for  the 
quadratic  given  by  equation  (2-53)   to  be  positive  definite. 
The  quadratic  will  assume  positive  values  for  any  variation 
in  P  if  and  only  if  the  determinant  of  D  given  by  equation 
(2-55)  and  all  its  principal  minors  are  greater  than  zero. 
For  symmetric  matrices,  the  first  positive  definite  condi- 
tion to  break  is  the  determinant  condition.     Hence,  the 
spinodal  is  defined  as  that  locus  of  points  satisfying 
det  D  =  0 . 

For  the  pure  component  case,  the  matrix  in  equation 
(2-55)  has  only  one  element  and  the  stability  criterion 
derived  from  it  is 

>  0  or  <  >  0  (2-56) 

T 

Equation   (2-56)   is  the  mechanical  stability  criterion. 

There  is  also  a  thermal  stability  criterion,  but  it  cannot 


3  2* 


3V 
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be  derived  using  the  Helmholtz  free  energy  representation. 
From  the  internal  energy  representation  we  find  that 


>  0 


or 


V 


(2-57) 


Equations   (2-56)   and  (2-57)  are  the  classical  stability 
criteria  whose  physical  content  is  known  as  the  Le  Chatelier 
principle.     According  to  this  principle,  the  criterion  for 
stability  is  that  all  spontaneous  processes  induced  by  any 
deviation  from  equilibrium  will  be  in  the  direction  of  re- 
storing the  system  to  equilibrium. 

For  the  case  of  a  binary  mixture  the  positive  definite 
condition  for  the  matrix  equation  (2-55)  is 


3V3N^ 

3N^ 

(2-58) 


with  all  its  principal  minors  greater  than  zero. 

Now,  it  is  possible  to  link  the  condition  given  by 
equation  (2-58)  with  the  equations  given  by  the  fluctuation 
theory.     Using  equations   (2-23)   and  (2-26)  we  get  for  the 
elements  of  the  determinant  in  equation  (2-58) 


3P 

3V 


T,N 
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1 
V  < 


^    [x^(l-C^^)   +  2x^X2(1-0^2 


+  X2(l-C22)] 


(2-59) 


9V 


) 


T,V,N2 


T,N 


3P_ 

8N, 


T,V,N2 


=     ^  [x^(l-C^^)    +  X2(l-C^2)^ 


(2-60) 


9N? 


T,V,N2 


3U] 


RT     r  _1_ 

V   ^p^  "  p 


(2-61) 


For  symmetric  matrices,  the  first  positive  definite 
condition  to  break  is  the  determinant  condition.  Hence, 
the  spinodal  for  a  binary  system  is  defined  as  that  locus 
of  points  satisfying 


3^A 


2 

3N7 


=  0 


(2-62) 
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Using  equations  (2-59) ,    (2-60) ,  and  (2-61)   in  equation 
(2-62) ,  we  get  the  following  relation  between  the  direct 
correlation  function  integrals  at  the  spinodal  curve 


^1^11  ■'"^2^22  "^1^2  ^^11^22  "^12^  -1 


(2-63) 


For  a  binary  system,  the  stability  criterion  when 
working  with  pressure,  temperature,  and  mole  fraction  as 
independent  variables  is  (Prausnitz,  1969) 


3xJ 


>  0 


(2-64) 


TP 


where  g  is  the  molar  free  energy  of  Gibbs.     Equation  (2-64) 
also  leads  to  equation  (2-63) .     A  similar  expression  is 
obtained  for  multicomponent  systems. 


3x^ 


^2 

_3_2_ 


3x     T  3x, 
n-1  1 


3x^  i3x, 
n-i  1 


^2 


3x 


n-1 


(2-65) 


with  all  its  principal  minors  greater  than  zero. 

It  might  be  simpler  to  check  stability  with  the  Gibbs 
free  energy,  but  then  (n-1) ^/2  +  (n-l)/2  second  derivatives 
must  be  evaluated  to  check  the  main  determinant.  With  the 
Helmholtz  free  energy,  for  a  given  pressure,  temperature 
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and  composition,  the  framework  of  the  fluctuation  theory 
can  be  used  directly.     The  only  difficulty  is  that  the 
density  of  the  solution  must  be  found  using  equation  (2-29)  . 

For  the  two  pure  component  points  in  the  mixture, 
equation  (2-63)   reduces  to 

C^i    =     1  (2-66) 

and  with  the  aid  of  equation  (2-59)  we  can  say  that  in  a 
pure  component  the  spinodal  is  defined  as  that  locus  of 
points  satisfying  the  limit  of  mechanical  stability. 

In  practice,  equation  (2-62)  may  be  of  little  utility 
because  equations  for  multicomponent  mixtures  have  many 
more  terms  to  handle. 


CHAPTER  3 

DIRECT  CORRELATION  FUNCTION  INTEGRALS  FROM 
EXPERIMENTS  AND  CORRELATION 


3 . 1  Introduction 

A  major  motivation  for  studying  the  fluctuation  solu- 
tion theory  approach  is  that  direct  correlation  function 
integrals  may  be  relatively  simple  and  universal  functions 
of  number  density  (or  density  at  macroscopic  level) .  A 
general  microscopic  approach  to  prove  this  point  is  not 
our  intention.     Rather,  we  examine  macroscopic  experimental 
information  about  the  integrals.     It  is  worth  noting  that 
the  relation    of   the  integrals  to  experimental  information 
is  through  derivatives  of  thermodynamic  properties.  As 
will  be  seen  in  Section  3.2  the  integrals  are  related  to 
the  isothermal  compressibility,  partial  molar  volumes  of 
components  in  a  mixture,  and  second  derivative  of  the 
excess  Gibbs  free  energy.     This    kind   of  dependence  can 
make  evaluation  somewhat  difficult. 

Section  3.2  presents  equations  and  results  for  ideal 
mixtures.     Section  3.3  does  the  same  with  nonideal  mixtures. 
Finally,  Section  3.4  shows  how  changing  parameter  values 
in  excess  Gibbs  free  energy  models  affect    the   form  of 
the  direct  correlation  function  integrals. 
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3 . 2    Experimental  Behavior  of  Ideal  Mixtures 

Thermodynamic  properties  of  binary  systems  are  com- 
pletely defined  by  C^^^,  C^2'  ^22'     Equations  (2-32)  and 
(2-26)  are  the  starting  point  through  the  expressions 


3  In  Yi 


(1  -  C.  .) 


i,j  =  1,2  (3-1) 


8  P/RT 


1  -     Z     x.C. . 
i=l  ^ 


(3-2) 


We  want  to  find  an  expression  for  the  composition 
variation  of  the  direct  correlation  function  integral. 
Operating  on  equations  (3-1)  and  (3-2)  we  get 


V  2  2 

icRT  -  ^  "  ^1^11  ~  ^2^22  "  ^^1^2^12 


(3-3) 


KRT  ^1^11  "  ^2^12 


(3-4) 


3  In  Y- 
S  x,' 


T,P 


=  X2  f(l-C^^)  (I-C22)  -  (1-C^2)^^^  (3-5) 


We  can  rearrange  these  last  equations  to  obtain 


V?  d  In  Y. 

(1  -  C     )   =   =—  +  X   =• 

"-11^       vkRT       ^2  dx. 


(3-6) 


T,P 
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(1-C3_2) 


(I-C22) 


VkRT  "  ^1 


^2 

VkRT  ^1 


d  In  Y]_ 


dx. 


(3-7) 


T,P 


d  In  Y2 


dx. 


T,P 


(3-8) 


The  experimental  behavior  of  the  direct  correlation 
integrals  can  be  determined  through  simultaneous  knowledge 
of  isothermal  compressibility,  molar  volume,  partial  molar 
volume,  and  partial  derivative  of  the  logarithm  of  the 
activity  coefficient.     This  is  the  pure  component  molar 
voliame  as  a  function  of  pressure  and  the  excess  Gibbs  free 
energy  as  a  function  of  composition  and  pressure.  The 
excess  Gibbs  free  energy  is  related  to  the  activity  coef- 
ficient by 


E 


,E 


^  RT 


=    X,   In  Yt  +  x-j  In  Y' 


(3-9) 


where  the  g  values  are  generally  found  from  a  model  whose 
parameters  are  obtained  from  phase  equilibrium  and  density 
data.  Taking  the  derivative  with  respect  to  mole  fraction 
of  component  1  twice,  we  get 


8x? 


T,P 


1  3  In  Yj^ 
^2  ^^1 

1    3  1^  ^2 


X, 


3x, 


T,P 


T,P 


(3-10) 


(3-11) 
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The  partial  molar  volumes  are  obtained  from  the  excess 
volume  by 


Vl  =  X,  II 


2         2         1  3x 


T,P 


T,P 


(3-12) 


(3-13) 


The  isothermal  compressibility  is  related  to  the  excess 
molar  volume  by 


1  8v^ 


V  3P 


=  < 


(3-14) 


T,x, 


where  <  is  related  to  the  solution  isothermal  compressibil- 
ity by 


o^E 

<    =  <,    +  +  K 

1  ^2  ^ 


(3-15) 


with  $^    as  the  molar  volume  fraction. 

For  an  ideal  solution  equations  (3.6)  to  (3.8)  become 
2 


(1-ct?) 


(1-ci^) 


(1-ci?) 


V, 


11  v<RT 


^2  2 
„o  o 

^1  ^2 

VKRT 

2 

o 

V-, 


22'  vkRT 


(3-16) 


(3-17) 


(3-18) 
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These  equations  show  that  at  constant  pressure  and 
constant  temperature  the  inverse  of  the  direct  correlation 
function  integrals  is  a  linear  function  in  mole  fraction. 

Table  3-1  presents  the  direct  correlation  function 
integrals  for  n-heptane — n-octane,  chosen  as  representative 
of  ideal  mixture,  at  different  temperatures.     A  good  indica- 
tion of  solution  ideality  is  the  agreement  between 
experimental  and  calculated  using  equation   (3-19) . 

Since  this  indicates  correlation  within  experimental  error, 
the  n-heptane — n-octane  system  is  ideal.     Figure  3-1  shows 
the  shape  of  the  direct  correlation  function  integrals 
which  are  not  linear  in  mole  fraction.     Figure  3-2  confirms 
the  ideality  of  the  system  showing  a  linear  variation  of  the 

inverse  of  (1-C. .)  with  mole  fraction.     Table  3-2  shows 

ID 

that  the  same  low  pressure  behavior  is  observed  at  high 
pressure.     Tables  3-3  to  3-10  present  direct  correlation 
function  integrals  for  different  ideal  systems.     In  all 
these  data,  differences  between  experimental  and  C22 

calculated  are  probably  within  experimental  error   (which  is 
0.5-3.%  in  <) .     Thus  we  can  conclude  that  these  systems 
are  essentially  ideal.     Another  possible  indicator  of 
ideality  would  be  the  correlation  coefficient  for  linear 
fitting,  r,  of  — Y        ^  function  of  x^.     In  the  above 
systems,  r  is  very  close  to  unity,  containing  ideal  behavior. 
However,  as  shown  below,  such  correlations  also  arise  in 
nonideal  systems, -meaning  that  linear  correlation  is  not 
useful  in  discriminating  nonideality. 
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igure  3-1.  Direct  correlation  function  integrals 
(DCFI)  in  n-heptane  (1) — n-octane  (2) 
at  298.15  K  and  1  atm  (•,  x,  +  data 

from  references  in  Table  3-1,   

smoothed) . 
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3.3    Experimental  Behavior  of  Nonideal  Mixtures 

We  have  chosen  the  acetone-carbon  disulphide  system  as 
representative  of  a  nonideal  mixture. 

Winnick  and  Powers  (1966a, b)  made  a  detailed  study  of 
this  mixture  at  low  and  high  pressure.     There  is  informa- 
tion to  calculate  excess  Gibbs  free  energy,  excess  molar 
volxime,  and  excess  isothermal  compressibility.     With  these 
data  and  equations  (3-6)   to  (3-8)  we  can  obtain  the  direct 
correlation  function  integrals  for  this  system.     Table  3-11 
presents  results  at  atmospheric  pressure  and  a  temperature 
of  273.15  K  using  Winnick  and  Powers'  equations.     For  the 
free  energy  at  atmospheric  pressure,  they  used  the  vapor 
pressure  data  of  Zawidsky  (1900)  and  a  fourth  order  equation 
to  smooth  the  free  energy  values.     In  order  to  convert  the 
data  from  298.23  K  to  273.15  K,  they  used  the  calormetric 
data  of  Schmidt  (1926)   and  Staveley  et  al.    (1955).  The 
correlation  coefficient  of  the  linear  fitting  of  the  in- 
verse of  (l"C-]^-j^)   is  almost  one  indicating  behavior  near 
ideal.     But  the  difference  between  the  experimental  and  the 
calculated  direct  correlation  function  integrals  for  a  pair 
of  carbon  disulphide  molecules  is  large  (from  20  to  100%) . 
Figure  3-3  shows  the  particular  shape  of  C^^,  C^2'  ^22- 
Figure  3-4  presents  the  inverse  of  (l-C^^)  which  shows 
that  the  linear  fitting  is  not  really  good.     The  points  do 
not  form  straight  lines  at  all.     The  extreme  points  of  the 
linear  fitting  at  mole  fraction  zero  or  one  are  not 
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Figure  3-3.     DCFI  in  acetone   (1) -carbon  disulphide  (2) 
at  273.15  K  and  1  atm  using  equations  of 
Winnick  and  Powers   ( 1966a, b)  . 
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Figure  3-4.     Inverse  DCFI  in  acetone   (1) -carbon 

disulphide   (2)   at  273.15  K  and  1  atm 
using  equations  of  Winnick  and  Powers 
(1966a, b). 
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coincident  with  the  values  of  pure  component  calculations. 
For  example,  at  acetone  mole  fraction  equal  to  zero,  the 
experimental  value  is  2.98E-2  and  the  linear  regression 
value  is  2.71E-2.     The  conclusion,  already  mentioned,  is 
that  the  correlation  coefficient  alone  is  not  a  good  indi- 
cator of  nonideality.     On  the  other  hand,  the  differences 
between  the  experimental  and  calculated  indicate  that 

the  mixture  is  not  ideal. 

It  is  apparent  from  the  curve  in  Figure  3-3  that  direct 
correlation  function  integrals  are  not  likely  to  be  modelled. 
It  is  possible  that  a  problem  could  be  in  the  experimental 
isothermal  compressibility.     Figure  3-5  shows  the  smoothed 
curve  of  Winnick  and  Powers  (1966a)   for  the  isothermal  com- 
pressibility as  function  of  mole  fraction  used  to  evaluate 
j .     There  is  a  discrepancy  at  the  pure  acetone  side  that 
could  change  the  form  of  the  curve.     An  alternative  fitting 
of  the  data  would  change  the  results  near  the  pure  acetone 
side.     However,  any  such  change  would  be  very  small  near  the 
pure  carbon  disulphide  side  where  C^^^^  goes  through  an  un- 
expected maximum. 

The  above  values  were  obtained  using  the  polynomial 
fit  of  Winnick  and  Powers   (1966b)   to  experimental  Gibbs 
free  energy.     Another  way  to  do  this  is  with  a  model. 
Since  the  acetone-carbon  disulphide  system  is  close  to  but 
not  at  immiscibility ,  the  Wilson  equation  is  appropriate 
to  describe  the  nonideality,  and  the  results  of  this  method 
will  be  adopted  as  representative  of  this  system.  Table 
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Figure  3-5.     Experimental  isothermal  compressibilities 
of  acetone   (1) -carbon  disulphide   (2)  at 
273.15  K  and  1  atm. 
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3-12  presents  values  of  derivatives  of  activity  coeffi- 
cients using  Winnick  and  Powers'  polynomial  and  the  Wilson 
equation.     Table  3-13  and  Figure  3-6  present  results  for 
the  same  calculation.     There  is  no  maximum  in  C^^,  and 
the  shape  resembles  that  of  ideal  mixtures.     Figure  3-7 
presents  the  inverses  of  the  (l-C^^).     The  high  values  of 
the  derivatives  at  infinite  dilution  explain  the  differ- 
ence in  C^j  obtained  using  equations   (3-6)   to  (3-8)  with 
Winnick  and  Powers'  data.     Before  discussing  the  conse- 
quences of  this  for  further  modelling,  we  show  the  results 
of  many  more  nonideal  systems  in  Tables  3-14  to  3-29. 

These  results  show  two  important  features.  First, 
there  are  some  systems  with  extrema  in  the  C^^ .  Second, 
there  are  some  systems  where  the  C^^  values  cross  as  a 
function  of  composition. 

3.4    Nonideal  Behavior  from  Excess  Gibbs 
Energy  Models 

As  shown  above,  the  results  for  the  C^^  are  sensitive 
to  the  model  used  to  describe  nonideal  behavior  of  mixtures 
This  section  gives  some  analytic  results  for  common  g^ 
models. 

The  simplest  model  for  g    is  a  two-suffix  Margules 
equation;  mixtures  that  follow  this  equation  are  known  as 
quadratic  mixtures.     The  starting  equation  for  binary 
mixtures  is 
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Figure  3-6.     DCFI  in  acetone   (1) -carbon  disulphide 
(2)   at  273.15  K  and  1  atm  using  the 
Wilson  Equation  for  excess  Gibbs 
energy . 


61 


5.0 


4.0 


3.0 


^22 

C^2  from  Figure  3-6 
^11 

Linear  regression 


2.0 


1.0 
0.0 

Figure  3-7 


JL 


0.25 


± 


0  .50 


0  .75 


1.00 


Inverse  DCFI  in  acetone   (1) -carbon 
disulphide   (2)   at  273.15  K  and  1  atm 
using  Wilson  Equation  for  excess  Gibbs 
energy . 
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^  RT 


-         X2  A 


(3-24) 


From  equation  (3-24)  and  using  thermodynamic  relations 
we  can  get  expressions  for  excess  entropy,  excess  enthalpy, 
and  excess  volume  as  given  by 


s  = 


-X      X  ^ 

^1  ^2  aT 


(3-25) 


h^  = 


^1  ^2  81/T 


(3-26) 


E  3A 
^    =  ^1  ^2  9P 


(3-27) 


Assximing  some  form  for  A  such  as 


A  =  A     +  A  P  +  A^T 
op  t 


(3-28) 


where  A^,  A^,  and  A^  are  constants,  we  get  for  thermodynamic 
properties  the  following: 


(3-29) 


V      =   X-    Xo  A 
1     2  p 


(3-30) 


h     =  X,   x„  [A^  +  A  P] 
±    z      o  p 


(3-31) 
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=  X,        A    +        P  (3-32) 


1     2  o 


=  X,  Xt  A  (3-33) 

1      Z  O 


KV  =  <^    +   ^2  ^2  ^2  (3-34) 


V,  =  v°      Ap  x2  (3-35) 


In  Y-  =  A  X?  (3-36) 


1  D 
d  In  y. 

— ^   =  -2A  X.  (3-37) 

d  x^  : 

In  order  to  simplify  further  our  equations,  we  assume 
that  A  is  a  constant  independent  of  pressure  and  temperature; 
then  for  a  regular  quadratic  mixture  (regular  because  ex- 
cess entropy  is  zero,  quadratic  because  Gibbs  energy  is 
quadratic  in  mole  fraction) ,  the  direct  correlation  func- 
tion integrals  are 

^o2 

-^11^  =   —  -  2^  4  ^^-'^^ 


■'^12^   =  —  ^  2A  x^  X2  (3-39) 


RT(Xj^  v°  K°  + 

X2 

o 
^2 

K2) 

O  O 

^1  ^2 

RT(Xj^  v°  <°  + 

X2 

o 
^2 

K2) 

o2 

^2 

RT(Xj^  v°  K°  + 

X2 

o 
^2 

K2) 

(1-^22)   =  —  -  2A  x2  (3-40) 
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In  the  particular  case  of  infinite  dilution  equations 
(3-38)  to.  (3-40)  reduce  to 

2 

<»1 

(l-^ll)       =  o    o    -  (3-41) 

RT  V2  K2 

o 

(1-C,2)"'  =  (3-42) 
RT  <° 

(1-C„)°°1  =   V  (3-43) 

22  RT  4 


where  the  superscript  "l  indicates  that  mole  fraction  of 
component  one  goes  to  zero. 


o 

(1-C,,)"2  =  ^  (3-44) 
RT 

(1  -CtO"'^  =    (3-45) 

RT  4 

2 

o 

(1-C,,)"^  =   ^— —  -  2A  (3-46) 

22  RT  v°  <° 


where  the  superscript  oo2  indicates  that  mole  fraction  of 
component  two  goes  to  zero. 

To  test  if  a  solution  modeled  by  the  Margules  equation 
could  show  maxima  and  crossing  of  the  C^^ ,  we  apply  these 
equations  to  the  acetone-carbon  disulphide  system.  A 
crossing  of  the  C^^  exists  only  if  the  unrealistic  value 
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of  4.45  is  used  for  A  (a  value  of  1  for  A  indicates  the 
limit  of  inmiscibility)  .     Further,  no  maximxam  can  be  found 
for  C^^  at  zero  acetone  mole  fraction,  unless  A  is  of  the 
order  of  5.00  E13. 

Another  empirical  equation  for  excess  Gibbs  free 
energy  (Abott  and  Van  Ness,  1975)   is  given  by 


x^  X2      ^21  ^1      ^12  ^2 


(A22^  ^12  ^1  "2 


(3-47) 


This  equation  with  its  four  constants  related  to  in- 
finite dilution  behavior  allows  more  variation  in  the 
behavior  of  direct  correlation  function  integrals.  The 
four  parameters  in  equation  (3-4  7)  are  given  by 


^12  =  1^  ^ 


ool 


(3-48) 


^12  ~  ^21  "  ^^12 


'*'21  "  ^^12  ^^21 


1  d  In  Yi 

2  dx^ 

1  d  In  Y2 

2  6kZ 


oo2 


(3-49) 


(3-50) 


(3-51) 


Table  3-3  0  presents  values  for  the  partial  derivative 
of  activity  coefficient  for  the  acetone-carbon  disulphide 
system  with  various  sets  of  infinite  dilution  values  in 
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equations  (3-4)   to  (3-51)  .     Figure  3-8  shows  C^^^^  values 
over  the  entire  concentration  range  for  the  original 
Winnick  and  Powers  values  and  for  the  arbitrarily  chosen 
infinite  dilution  derivatives.     Figure  3-9  shows  the  dilute 
acetone  (x^  -»■  0)  behavior  of  C^^^  for  all  values  used  in 
Table  3-30.     Finally,  Figure  3-10  presents  the  shape  of 
the  excess  Gibbs  free  energy  for  each  set  given  in  Table 
3-30. 

In  conclusion  we  can  say  that  simple  forms  of  C^^ 
reproduce  very  well  the  excess  Gibbs  free  energy  calculated 
with  complicated  forms  of      ^  .     In  other  words,  the  figures 
indicate  that  modelling  the  C^^  would  be  desirable  for  pre- 
dieting  g    values  because  results  are  not  extremely  sen- 
sitive to  C^j  values,  particularly  at  high  dilutions.  On 
the  other  hand,  the  same  insensitivity  indicates  that  the 
"experimental"  direct  correlation  function  integrals  are 
uncertain.     This  uncertainty  in  the  C^^  is  a  major  obstacle 
to  modelling. 
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Figure  3-8.     -Cin   for  acetone   (1)   in  carbon  disul- 

phxde  at  273.15  K  and  1  atm  from  various 
expressions  for  G^ . 
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Figure  3-9.     -C]^]_  for  acetone   (1)   in  carbon  disulphide 
at  27  3.15  K  and  1  atm  from  various  ex- 
pressions for  . 


Figure  3-10.  Excess  Gibbs  energies  for  acetone  (1) 
and  carbon  disulphide  at  273.15  K  and 
1  atm  from  various  expressions. 


CHAPTER  4 

MODELS  FOR  DIRECT  CORRELATION  FUNCTION  INTEGRALS 

4 . 1  Introduction 

In  order  to  be  applied,  the  fluctuation  solution  theory 
presented  in  Chapter  2  requires  a  model  for  the  direct  cor- 
relation function  integrals.     This  chapter  describes  ex- 
perience with  one  such  model.     Its  concept  is  that  condensed 
fluid  behavior  can  be  expressed  as  that  of  a  hard  sphere 
fluid  plus  a  perturbation.     Here,  the  equation  for  the 
direct  correlation  function  integrals  is  a  hard  sphere 
expression  plus  a  linear  density  term.     Section  4.2  de- 
scribes this  idea  as  used  by  Mathias  (1978)   for  gas  solu- 
bility.    Section  4.3  introduces  a  modification  in  the  hard 
sphere  expression  to  more  closely  describe  pure  liquid 
compressions.     Section  4.4  investigates  the  possibility  of 
replacing  the  hard  sphere  equation  with  a  hard  convex  ex- 
pression.    Finally,  Section  4.5  shows  application  of  the 
Mathias  model  to  solubilities  of  gases  in  binary  solvents. 

4 . 2    Description  of  the  Basic  Model 

The  basis  of  the  present  model  for  direct  correlation 
function  integrals  comes  from  work  by  Mathias  (1978,  1981). 
We  will  describe  its  physical  picture,  generalize  its 
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equations,  present  its  extension  to  multicomponent  systems, 
and  explain  its  working  procedure. 

In  the  early  1970s,  it  was  shown  that  the  structure  of 
a  simple  pure  fluid  is  basically  the  structure  of  an  equivalent 
hard  sphere  fluid  (Chandler  et  al.,  1983).  Intermolecular 
forces  are  separated  into  short-range  repulsion  and  long- 
range  attraction.     At  high  densities  the  structure  of 
liquids  is  determined  by  packing  effects  related  to  the 
short-range  repulsion.     On  the  other  hand,  the  long-range 
attraction  generates  a  uniform  background  potential  that 
keeps  the  liquids  at  its  high  density  but  does  not  affect 
the  structure. 

For  hard  spheres,  analytical  expressions  such  as  the 
Percus-Yevick  compressibility  equation  (Lebowitz  and 
Rowlinson,  1964)  give  accurate  compressibility  factors  at 
densities  up  to  three-fourths  of  that  of  the  solid-liquid 
transition.     For  higher  densities,  the  semiempirical  ex- 
pressions such  as  Carnahan-Starling  (1969)  equation  are 
better.     As  a  result,  the  thermodynamic  properties  of  hard- 
sphere  fluids  are  well  described  over  the  full  range  of 
density. 

The  properties  of  real  spherical  liquids  may  be  calcu- 
lated by  a  perturbation  technique  with  the  hard  sphere  as 
a  reference  and  with  the  attractive  forces  as  a  perturba- 
tion.    For  the  direct  correlation  function  integral  of  a 
pure  liquid  the  simplest  expression  of  this  physical 
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picture  means  a  hard  sphere  expression  plus  a  linear  density 
correction  term 


where  p  is  the  density  reduced  by  a  suitable  spherical 
covolxame. 

Before  giving  particular  functional  forms  and  variables 
in  equation  (4-1) ,  it  should  be  noted  that  Gubbins  and 
O'Connell   (1974)   and  Brelvi  and  O'Connell   (1972)   showed  that 
a  corresponding    states   behavior  can  be  found  for  the 
direct  correlation  function  integrals  of  pure  components. 
Figure  4-1  shows  the  generalized  correlation  for  C  of 
Brelvi  and  O'Connell   (1972) .     The  correlation  works  quite 
well  for  densities  larger  than  1.5  p*,  where  p*  is  a  charac- 
teristic parameter  close  to  the  critical  density  for  non- 
polar  substances  but  larger  for  polar  species.  Mathias 
(1978)   introduced  a  two  parameter  corresponding  states 
correlation  in  the  form  of  equation  (4-1) . 


C(T,p)   =  C     ip)   +  p  B(T) 


(4-1) 


C(T,p)   =  C^^(p,T)    -  2   p[g(T)    -  f(T)] 


(4-2) 


with 


C^^(p,T)    =  $(n)    =  $[|  f (T) ] 


(4-3^) 


where  the  hat  "  indicates  reduced  property,  and  n  is  the 
packing  fraction,  here  a  function  of  reduced  temperature 
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Generalized  correlation  for 
pure  component  DCFI  with  re- 
duced density   (Brelvi  and 
O'Connell,   1972) . 


91 


and  density.     Equation  (4-2)  uses  two  parameters,  a  charac- 
teristic temperature,  T*,  and  a  characteristic  volume,  V*, 
or  density,  p*.     The  two  functions  f  and  g  were  determined 
by  Mathias  from  data  on  pure  argon.     Appendix  A  presents  a 
complete  list  of  his  equations  and  functions. 

Introducing  equation   (4-2)   in  equation   (2-29)  we  ob- 
tain the  following  equation  of  state  for  a  pure  liquid  phase 

=  P^d  +        -  2nf)/(l  -  rif)^ 

-  P^'d  +        -  2n^)/(l  -  rij-)^ 

+   (pf  -  pj)   V*    (g(T/T*)   -  f(T/T*))  (4-4) 

where 

T]  =  ^  f  (T/T*)  (4-5) 

The  correlation  of  pure  compression  data  is  generally 
good,  with  the  results  being  sensitive  to  the  value  of  V* 
but  insensitive  to  T*.     As  our  objective  is  to  work  with 
mixtures,  the  main  use  of  equation   (4-4)   is  only  to  obtain 
pure  component  parameters. 

For  mixtures  it  is  assumed  that  the  pair  direct  cor- 
relation function  integrals  are  given  by  an  extension  of 
equation  (4-2) 
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C.  .  (T,d)  =  C^?  -  2pV* .   [g.  .  -  f .  .]  (4-6) 
ID     ^  ID  ID       ID  ID 


with 


C^^.  =  $(£)  (4-7) 


where 


TT       ^  .  V*.f. 

^      6  Pi   ^2  "n     )     k  =  0,  1,  2,  3  (4-8) 


and 


f.  =  f(T/Tt.)  (4-9a) 
1  11 

in  addition 

g^j  =  g(T/Tjj)  (4-9b) 

and 

jl/3  ^  ^1/3  2 
fij  =   (   )  (4-9C) 

In  previous  work  (Mathias,  1978) ,  the  hard  sphere 
mixture  direct  correlation  function  integral,  cV?,  is 
described  by  the  Carnahan-Starling  equation  for  mixtures. 
Mixing  rules  are  required  for  evaluating  the  characteris- 
tic parameters,  V^^  and  Tt^  for  i^^ j  .     The  mixing  rules 
used  are 
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where  V^^  is  the  characteristic  volume  for  pure  i 


1/2 

T*  .   =    (T*.   T*  .)  (1  -k.  .)  (4-11) 


where  T*^  is  the  characteristic  temperature  for  pure  i,  and 

is  the  binary  parameter  which  must  be  evaluated  from 
binary  data. 

Introducing  equation  (4-6)  in  equation  (2-29) ,  we 
obtain  the  following  equation  of  state  for  a  multicom- 
ponent  liquid  phase 


-  ^     P^^T.pf)    _  P^^(T,  £^) 

RT  RT  RT 


n      n         _  _ 
+     11     [p.p     -  P^P^]  V*.  (g..-f..) 

j=i  ^i^D     ^i^:-*    13  ^iD  1: 


(4-12) 

Similarly,  introducing  equation  (4-6)   in  equation  (2-3  0) 
we  can  obtain  an  expression  for  activity  coefficients  of 
multicomponent  liquid  mixtures  (see  Appendix  A  for  details) . 

4 . 3    Comparison  of  Different  Hard  Sphere  Equations 

Mathias  (1978)  used  the  Carnahan-Starling  equation  for 
the  reference  term  because  it  is  the  best  analytic  expres- 
sion for  hard  spheres.     This  fixed  the  reference  term  in 
equation  (4-1) ,  and  he  found  a  functional  fom  for  the  hard 
sphere  diameter  and  the  linear  densitv  coefficient  from 
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compression  data  on  argon,  the  reference  compound.     An  al- 
ternative to  his  procedure  is  to  keep  free  the  reference 
equation.     Telotte  and  O'Connell   (1982)  used  this  approach 
in  fitting  pure  component  data  to  a  generalized  hard  sphere 
equation  of  state 

^hs  ^  1  +  ^  +        +  an^  (4.13) 
(1  -n) 

or  its  pressure  equivalent 

(1  -  n) 


or 


,hs 


(p,T,a)  =  <I>(n/  a) 


(4-3b) 


where  a  is  a  parameter,  and  it  is  related  with  the  steep- 
ness of  isotherms  in  the  P-p  diagram.     Different  choices 
of  a  yield  different  hard  sphere  equations.     The  Carnahan- 
Starling  expression  results  from  a  value  of  a  =  -1,  the 
Percus-Yevick  compressibility  equation  from  a  values  of 
a  =  -3,  and  the  Percus-Yevick  virial  equation  from  a  =  0. 

Telotte 's  approach  was  to  treat  a  as  adjustable 
parameter  in  the  building  of  the  model.     He  chose  methane 
as  the  reference.     Minimizing  the  difference  between  the 
experimental  direct  correlation  function  integral  given 
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by  equation  (4-2)  with  C  from  equation  (4-3b)  ,  an  optimxim 
value  of  -4.2  was  found  for  a. 

The  choice  of  hard  sphere  equations  to  use  as  a  ref- 
erence is  critical  not  only  for  pure  components,  but  also 
in  mixtures.     Equations  with  the  adjustable  parameter  a  for 
mixture  pressure  and  direct  correlation  function  integrals 
in  a  multicomponent  mixture  are  presented  in  Appendix  B. 
We  describe  our  analysis  with  the  n-heptane — n-octane 
system  (ideal  mixture)  and  with  acetone-carbon  disulphide 
system  (nonideal  mixture)  to  determine  the  sensitivity  of 
the  composition  dependence.     A  simple  way  to  do  this  is  to 
evaluate  the  hard  sphere  diameter  for  each  component  when 
pure.     At  fixed  a  the  diameter  is  calculated  by  taking  the 
hard  sphere  direct  correlation  function  integral  equal  to 
the  corresponding  experimental  value.     Then  the  differences 
between  the  experimental  data  and  the  calculated  values 
using  equation  (B-6)  are  examined  over  the  entire  concen- 
tration range.     Table  4-1  presents  results  for  n-heptane 
with  n-octane  at  constant  pressure.     The  B^^  column  is 
the  difference  between  the  experimental  and  the  calculated 
C^j  divided  by  the  solution  density.     Figure  4-2  displays 
the  experimental  and  the  calculated  direct  correlation 
function  integrals  for  the  pairs  of  n-heptane.     The  small 
influence  of  varying  the  a  value  is  clear  from  the  figure. 
The  same  is  true  for  the  other  direct  correlation  function 
integrals  in  Table  4-1.     Further,  while  the  value  of  B. . 
varies  strongly  with  composition  in  all  cases,     the  value 
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0  .0  0.2  0  .4  0.6  0.8  1.0 

Figure  4-2.     Isobaric   (1  atm)   -Cn  values  in  n-heptane 
(1) — n-octane    at    333.15  K. 
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does  not.     Therefore,  the  compositional  variation 


in  mixtures  must  be  correlated  differently  than  these  of 
pure  fluids. 

Table  4-2  presents  results  for  n-heptane — n-octane  at 
a  constant  density  and  Figure  4-3  shows  direct  correlation 
function  integrals  for  the  n-heptane — n-heptane  pair. 
Again,   there  is  only  a  small  difference  between  the  num- 
bers from  the  different  hard  sphere  equations. 

Similar  results  for  the  system  acetone-carbon  disul- 
phide  are  presented  in  Tables  4-3,   4-4  and  Figures  4-4, 


Complex  molecules  have  nonspherical  hard  cores .  For 
modelling  of  direct  correlation  function  integrals  we  inves- 
tigate the  possibility  of  using  an  equation  for  a  convex 
body  of  revolution. 

The  direct  correlation  function  integral  for  hard 
convex  body  is  obtained  in  this  work  from  Boublik's  com- 
pressibility factor  equation  (Boublik,  1975) 


4-5  . 


4  . 4     Hard  Convex  Body  as  Reference 


^  _  1  +  (3y  -  2)  n  +  (3y^  -  3y  +  1)  - 
z  -  =  

(1  -  n) 


(4-15) 


where  z,  the  compressibility  factor  is 


z  = 


P 


(4-16) 


RTp 
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Figure  4-3.     Isochoric   (135.37  cc/mol)   -Cn  values 

in  n-heptane   (1) — n-octane  at  333.15  K. 
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 1  1 —   I 

Experiment   (Table  3- 

-13) 

—  —  —     From  equation  (B-6) 

- 

with  a  =  0,  -3 

- 

\  \ 

\  \ 

\  \ 

\  \ 

- 

\  \ 

V  \ 

\  \ 

\  \ 

\  \ 

\  \ 

\  \ 

\  \ 

\  X 

\  \ 

\  \ 

1                 1  1 
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0.0  0.2  0.4  0.6  0.8     x,  1.0 


Figure  4-4.     Isobaric   (1  atm)   -C^  values  in  acetone 
(1) -carbon  disulphide  at  273.15  K. 
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Figure  4-5.     Isochoric   (58.67  cc/mol)   -C^^  values  in 

acetone   (l)-carbon  disulphide  at  273.15 K. 
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and  Y  relates  the  hard  convex  body  voliime  (V)  ,  the  hard 
convex  body  surface  area  (S)  ,  and  the  (1/4  II) -multiple  of 
the  mean  curvature  of  the  hard  convex  body  (R)  by 


Y  =  If  (4-17) 


This  equation  shows  satisfactory  agreement  with  computer 
simulations  (Vieillard-Baron,  1974;  Rebertus  and  Sando, 
1977;  Nezbeda  and  Boublik,  1978;  Monson  and  Rigby,  1978) 
over  the  density  range  available.     The  Carnahan-Starling 
equation  for  hard  sphere  is  recovered  from  equation  (4-15) 
with  Y  =  1.     The  relation  between  the  direct  correlation 
function  integral  and  the  compressibility  factor  is  given 
by  the  following  expression 


1  -  c(p,T)   =  IfZM 

dp 


_  9pz 

9p 


(4-18) 


Using  equation  (4-15)  in  equation  (4-18) ,  and 


(1  -  C)   =   i — J  (1  +  (6Y-2)n 

(1-n) 


+  (9y^  -  6y  +  i)n^  -  4Yn^  + T^n^)  (4-19) 


Values  of  (1-C)   for  different  n  and  y  are  presented 
in  Table  4-5.     The  hard  convex  fluid  shows  less  compres- 
sibility with  larger  y •      As  complex  molecules  are  non- 
spherical,  the  use  of  a  hard  sphere  as  the  hard  core 
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reference  implies  an  adjustment  in  the  hard  sphere  voliame. 
The  corrected  volume  is  larger  than  the  real  volume  because 
the  compressibility  for  a  convex  body  is  smaller. 

Table  4-6  presents  the  comparison  of  this  adjustment. 
The  corrected  volume  was  chosen  to  minimize  the  difference 
between  the  (1-C)   for  a  hard  convex  body  and  the  (1-Cj  for 
a  hard  sphere.     The  adjustment  is  relatively  good  for 
n  <  .5  as  shown  in  Figure  4-6.     However,  at  higher  values 
of  the  agreement  is  poorer. 

Since  the  use  of  the  hard  convex  body  expression 
involves  another  parameter  which  probably  cannot  be  chosen 
from  molecular  structure,  given  the  relatively  satisfactory 
agreement  of  the  hard  sphere  equation  and  the  expectancy 
that  there  must  be  further  corrections  to  the  model  with 
perturbations,  the  hard  sphere  equation  appears  adequate 
at  this  stage. 

4 . 5    Applications  to  Gases  in  Mixed  Solvents 

Mathias  was  very  successful  in  using  the  model  de- 
scribed in  Section  4.2  in  applications  involving  liquids 
containing  supercritical  components.     This  extends  his 
analysis  to  calculate  Henry  constants  in  mixed  solvents. 

The  common  approach  to  calculate  Henry ' s  constants  in 
mixed  solvents  involves  Henry ' s  constants  in  the  pure 
solvents.     The  simplest  relation  was  proposed  by  Kritchevski 
(1937)  and  is  valid  for  an  ideal  behavior  of  the  liquid 
phase 
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Figure  4-6.     DCFI  for  pure  component  rigid 
bodies . 
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n 

In        ^  =     Z     X.   In         .  (4-21) 
±  /  in      j_=2  ^ 


where        ^  is  the  Henry  constant  of  component  1  in  a  mixture 

m,  and        .  is  the  Henry  constant  in  solvent  i. 
1/1 

If  the  liquid  phase  behaves  nonideally/  equation  (4-21) 
must  be  changed   (O'Connell,  1971b;  Prausnitz  and  Chueh/ 
1968) 


n 

In  Ht  ^  -    Z    X.  In  Ht  . 
'        i=2  ' 


n 

=  lim  Yi  „  -     I     X.    dim         .)  (4-22a) 
X  /in       1=2  1 ' 

x^  ->  0  -y  0 


or  in  unsymmetric  convention 


n 

In        _  -     Z     X.   In  H,  . 
Ifin  1  I'l 


n 


lim  V, 

'  1  ,m 

y  1 


Z  X.  dim  Yi  A ) 
i=2    ^  ^'^ 


(4-22b) 


X 


Depending  on  the  model  used  to  describe  the  activity  co- 
efficient different  forms  of  equation  (4-22a)  appear.  Two 
of  them  are  discussed  in  O'Connell's  work.     The  first  uses 
the  Margules  expression  and  transforms  equation  (4-22a) 
into 


Ill 


n  n  n-1 

In  H.  „  -     Z     X.   In  H-    .   =     I       I     a.,  x.x, 
l,in  ^  i>k  k=2  ^  ^ 


E 

~  ~  RT     (solvents  only)  (4-23) 


For  a  particular  case  of  a  solute  in  a  binary  solvent  equa- 
tion (4-23)   reduces  to 

In  H,  _  -  X-  In  Ht   „  -  x^  In        ,  =  -a^-,  x_x, 
l,m         Z  1,2  J  1,3  23     2  3 

E 

=  -  ^  (solvents  only) 
(4-24) 

The  second  derives  from  a  model  of  fluctuation  solution 
theory  presented  in  Chapter  2.     In  the  case  of  binary 
solvents  take  the  following  form 


In  -  Xo  In         _  -  x^  In  , 

l,m        z  1,2        J  1,3 


x^  9  In  Y2 
-T  9x3 


+  A-^  4)^(J)^  (4-25) 

T,P,N2 


where 


1050    |v?  -  v°l 

^23  =  ^1/2     o     o  (4-26) 
T  /     V2  V3 

4)^  is  the  voliame  fraction  of  solvent  i  in  pure  solvents 
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X.  v° 

(}).    =  i   (4-27) 

1  o  o 

^1^1  ^2^2 

Other  alternatives  are  to  change  the  mole  fraction 
in  equation  (4-22)   to  volume  fraction,  the  surface  fraction, 
or  by  using  the  r (volume)  and  q (surface)  parameters  of  the 
UNIQUAC  activity  coefficient  model   (Abrams  and  Prausnitz, 
1975)  . 

The  Mathias  method  uses  the  following  equation 


1"  «l,m  =  1^  «1,R  ^  1^  ^1  ('^'P'  =  Pm'   P""  =  Pr)  (^-28) 


where 


—  =     I    x^V  +  v^   (T,x^^)  (4-29) 
Pm      i=2    ^  ^ 

SF 

with  X      as  the  mole  fractions  free  of  solute. 

Equation  (4-28)   comes  from  equation  (A-4)   for  the  varia- 
tion of  solute  fugacity  from  infinite  dilution  in  one 
solvent  (R)   to  infinite  dilution  in  another  (mixture,  m, 
here) ;         is  calculated  using  equation   (A-4) .     The  binary 

constant      .  is  determined  from  k^^  and  data  on  H,    .  and 
li  IR  1,1 

to  insure  consistency  in  the  limit  x^      1 . 

Table  4-7  presents  results  comparing  Henry's  constant 
calculations  for  4  0  binary  solvent  systems.     The  numbers 
in  the  table  indicate  the  absolute  average  difference 


113 


0) 
H 
Xi 
(0 
Eh 


C 
•H 

to 

(U 
CO 
(0 

i-l 
o 

to 

4-) 

c 

(0 

tn 
C 
O 
u 

CO 

>i 

c 

0) 

C  m 
+J 

0)  (1) 
•P  > 

(0  iH 
-H  O 

s  w 

•H  >i 

(0  n 
U  (0 

c 

:  -H 


c 
o 

+J 

(0 

> 

Q 

(1) 

cn 

M 
0) 
> 
< 

H 

o 

(0 


o 
II 


00 

I 


CO 
I 


I 


tN 
I 

'3' 


CO 


(N 
I 


CM 
I 


I 


6 
d) 
JJ 
CO 
>i 
c« 


f*1 


O 


VD 
I 


'3' 

CO 


00 

I 


O  rH 


n 

CT> 

(N 

VD 

VD 

CTv 

CO 

VD 

VD 

(N 

O 

00 

(N 

i-H 

rH 

rH 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

rH 

O 

n 

rH 

iH 

r~ 

(N 

in 

00 

If) 

'3' 

(N 

CO 

in 

rH 

rH 

1 

+ 

+ 

+ 

+ 

+ 

rg 
cn 


00 


in  rH 


o 


in 


CO 


rH 

CO 

VD 

'3' 

CO 

CO 

in 

ro 

CO 

'3' 

rH 

1 

+ 

+ 

1 

rH 

VD 

CM 

in 

in 

rH 

c^ 

rH 

in 

rH 

CO 

ro 

VD 

ro 

ro 

+ 

00 

<N 

VO 

VO 

ro 

in 


ro 

CO 

(N 

ro 

VO 

ro 

PO 

o 

in 

CO 

rH 

CN 

in 

o 

rsi 

in 

rH 

O 

vD 

CN 

VO 

00 

PO 

c^ 

O 

vD 

00 

'3' 

00 

VD 

00 

VD 

rH 

ro 

(N 

ro 

rH 

ro 

rH 

rH 

rH 

H 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

1 

1 

+ 

ro 
VD 


ro 
+ 


"3- 

ro 

ro 

ro 

ro 

in 

CM 

ro 

VO 

CTi 

'3' 

cn 

o 

rH 

CTi 

CM 

00 

fN 

o 

o 

in 

CN 

ro 

ro 

in 

rH 

VD 

in 

VD 

rH 

ro 

rH 

rH 

rH 

ro 

00 

VD 

rH 

rH 

rH 

rH 

rH 

rH 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

1 

1 

+ 

VO 


CO 
VD 


cn 
ro 


in 


VD 

in 

vD 

CN 

CN 

CD 

rH 

VD 

00 

ro 

in 

VD 

CN 

rH 

in 

VO 

rH 

rH 

CO 

in 

rH 

rH 

rH 

CN 

ro 

CO 

CD 

rH 

O 

ro 

rH 

CN 

VD 

in 

rH 

CN 

rH 

rH 

rH 

rH 

rH 

rH 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

I 

1 

+ 

in 

00 

00 
'3' 

CM 

VD 

00 

cn 

rH 

CD 

in 
ro 

cn 

CO 

CM 

o 

in 

"3' 

r- 

CM 
CM 

rsi 

CN 
rH 

ro 

rH 

rH 

rH 
rH 

CO 
rH 

in 

rH 

00 

VO 

rH 
rH 

cn 

ro 

rH 

rH 

rH 

CO 

rH 

VO 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

1 

1 

1 

+ 

tn 

rg 

ro 

fN 

1 

in 

o 

'3' 

VO 

'3' 

CD 

CO 

CO 

rH 

VO 

rH 

VO 

O 

1^ 

rH 

rH 

CTi 

CD 

VD 

ro 

<N 

CN 

CN 

rH 

rH 

H 

rH 

rH 

ro 

rH 

VD 

00 

o 

CD 

n 

PO 

CM 

CD 

r» 

in 

in 

CD 

in 

rH 

VD 

in 

rH 

in 

r>- 

rH 

CN 

ro 

in 

CO 

ro 

rH 

CD 

CN 

o 

ro 

rH 

rH 

r- 

CN 

rH 

rH 

rH 

rH 

rH 

CN 

ro 

in 

VO 

00 

CD 

O 

rH 

PO 

rH 

rH 

rH 

rH 

H 

CN 

ro 


00 


CN 


VO 

H 

in 

CD 

VO 

VD 

PO 

• 

• 

• 

• 

• 

o 

o 

rH 

H 

rH 

CM 

CO 

CS) 

CN 

in 

O 

00 

CM 

in 

• 

• 

« 

• 

rH 

in 

o 

PO 

CN 

CN 

rH 

CN 

in 

in 

in 

PO 

o 

00 

CN 

CN 

ro 

00 

r« 

PO 

ro 

rH 

H 

vD 

CO 

CD 

o 

rH 

rH 

rH 

rH 

CN 

114 


00 
CM 


00 
I 


H 


o 

iH 

00 

r-l 

rH 

1-1 

ro 

rH 

cr> 

rH 

rH 

ro 

ro 

ro 

ro 

r^ 

rH 

VD 

VD 

ro 

O 

00 

vc 

o 

in 

00 

fN 

ro 

o 

O 

CO 

r~ 

CN 

rH 

o 

00 

•'T 

ro 

i-H 

<n 

CO 

CM 

rsi 

'a* 

O 

VD 

rH 

<Ti 

ro 

rH 

ro 

in 

vo 

1 

CN 
1 

1 

ro 
1 

ro 
1 

<N 
1 

CN 

1 

ro 
1 

CN 

CM 
1 

ro 

fN 

CN 

CN 

1 

CN 

fM 

fN 

ro 
1 

CM 
1 

ro 

00 

iH 

o 

( — 

ro 

CO 

rH 

fN 

rH 

rH 

in 

VO 

o 

00 

ro 

CN 

oo 

rM 

rH 

r- 

ro 

rH 

a\ 

rH 

ro 

ro 

in 

cn 

00 

CN 

1— 1 

CN 

CN 

ro 

r— 

in 

in 

in 

rH 

in 

in 

rH 

rH 

1 

+ 

+ 

+ 

+ 

+ 

+ 

00 

iH 

rH 

O 

\D 

ro 

fO 

rH 

ro 

ro 

cn 

C7^ 

O 

<-i 

H 

fN 

in 

00 

rM 

fN 

CM 

o 

CO 

rH 

CM 

o 

,—1 

in 

VD 

ro 

00 

in 

r~ 

rH 

in 

ro 

CTi 

CO 

o 

fN 

rH 

rH 

rH 

rM 

rH 

rH 

+ 

1 

1 

1 

CN 

ID 

ro 

CN 

ro 

0^ 

fN 

00 

in 

CN 

r- 

ro 

CO 

in 

ro 

rH 

O 

rs) 

<n 

ID 

"S- 

ro 

"a- 

CTi 

rH 

fN 

r- 

O 

>X) 

rH 

cn 

00 

VD 

H 

(Tv 

rM 

rH 

fN 

ro 

ro 

rH 

CO 

t-i 
1 

CN 

1 

rH 
1 

1 

ro 
1 

ro 
1 

ro 
1 

CN 
1 

CN 
1 

rH 

rH 

ro 
1 

CN 
1 

fN 

CN 
1 

rH 

CN 

fN 

1 

fN 
1 

CN 
I 

in 

in 

fN 

00 

VO 

rH 

r~ 

CM 

O 

00 

ro 

"a" 

"a" 

CN 

VD 

VD 

in 

o 

o 

in 

o 

o 

CO 

'a' 

ro 

cn 

O 

cn 

ro 

vc 

o 

ro 

o 

ro 

rH 

o 

CO 

\D 

C^ 

o 

in 

'a' 

lO 

m 

in 

in 

1 

(N 

1 

1 

rH 

1 

ro 
1 

CM 
1 

(N 

1 

fN 
1 

ro 
1 

rH 

rH 

fN 
1 

CN 

rH 

1 

rH 

rH 

rH 
1 

rM 
1 

rH 

1 

ro 

00 

ro 

00 

r~ 

00 

rH 

in 

in 

'a- 

rH 

rH 

in 

ro 

00 

o 

ro 

rM 

CM 

ro 

O 

o 

iH 

ro 

rH 

o 

00 

VD 

ro 

cn 

ro 

o 

00 

in 

vc 

in 

«£i 

r~ 

r- 

fN 

in 

00 

CO 

iH 

ro 

CN 

ro 

ro 

ro 

ro 

rM 

rM 

ro 

ro 

rM 

ro 

rM 

fN 

CM 

ro 

ro 

fN 
I 

'5' 


fN 
I 

"a- 


I 

'a" 


in 

VO 

o 

O 

CO 

VO 

ro 

T 

in 

ro 

00 

VO 

H 

rH 

CM 

rM 
1 

cn 

o 

00 

ro 

ro 

rH 

o 

rM 

in 

m 

CN 

rH 

rH 

r- 

in 

cn 

<n 

VD 

fN 

in 

VD 

o 

00 

VD 

VD 

rM 

rM 

o 

cn 

r~ 

o 

in 

m 

cn 

in 

O 

rH 

ro 

rH 

fN 

ro 

rH 

ro 

rH 

ro 

in 

rH 

CM 

ro 

rM 

VD 

CO 

cn 

rM 

ro 

r- 

ro 

'a' 

00 

rH 

00 

CM 

o 

cn 

VO 

"a- 

rM 

•a- 

rM 

cn 

rH 

rM 

rM 

in 

o 

00 

VD 

o 

in 

o 

o^ 

vD 

CM 

o 

CO 

o 

1 

rH 

rH 

ro 

VD 

ro 

rH 

VD 

rM 

in 

rM 

VD 

T 

VD 

VD 

rH 

in 

rH 

rH 

'a* 

rH 

rM 

rH 

rM 

rH 

.rH 

rM 

H 

rM 

rH 

rH 

rH 

rH 

rH 

rH 

in 

ro 

r- 

O 

in 

rH 

rM 

O 

ro 

CO 

ro 

rH 

in 

VD 

in 

H 

VD 

CM 
1 

H 

rH 

00 

rH 

CO 

in 

CO 

o 

VD 

O 

'a' 

o 

O 

r-~ 

en 

cn 

cn 

00 

CM 

'a' 

CN 

rH 

in 

rH 

ro 

r~- 

'a' 

VD 

VD 

CO 

rM 

O 

VD 

00 

CO 

cn 

in 

VD 

ro 

rH 

rH 

rH 

rH 

rH 

rH 

CM 

ro 

rH 

CM 

ro 

£ 

4J 

(0 

rH 

rM 

ro 

'a' 

in 

VD 

00 

cn 

O 

rH 

rM 

ro 

'a' 

in 

VD 

00 

cn 

O 

>l 

w 

rM 

CM 

rM 

CM 

rM 

r>) 

rM 

rM 

ro 

ro 

ro 

ro 

ro 

ro 

ro 

ro 

ro 

ro 

'3' 

115 


r3 
<U 
3 
C 
•H 
+J 

c 
o 
u 


I 

Q) 

rH 

(0 
E-" 


00 

I 


c 

•H 

§ 

m 

•rl 

4J 

>1 

0 

U 

H 

z 

(0 

+J 

f  r 

SOU] 

■p 

no 

u 

(0 

u 

0) 

x; 

0 

4-1 

(U 

to 

•  E 

tn 

c 

tn  n3 

0) 

0 

1  M 

4J 

•H 

^  (0 

(0 

4J 

u 

(0 

C 

•H 

•H 

• 

0  u 

> 

T3 

0 

c 

0) 

-0 

U 

U 

0) 

OJ  CO 

c 

0) 

-0 

> 

•H 

c 

■H 

•H 

CO 

C  (0 

(0 

4J 

c 

•H  JJ 

(0 

0 

10 

0) 

CT^ 

o 

(fl 

> 

<U 

C  h 

•H 

c 

c 

0 

4J 

0) 

■M  - 

•H 

>l 

0) 

■P  C 

CO 

iH 

ja 

(C  0 

0 

J-> 

3  -H 

ft 

to 

C  -U 

0 

0 

0)  U 

<: 

6 

c 

ro 

4J 

CO 

(0 

to 

c 

s: 

0 

•H 

•  CO 

116 


between  the  calculated  and  the  experimental  values.  Table 
4-8  contains  the  reference  of  the  data  and  the  systems 
involved . 

The  numbers  obtained  with  the  Mathias  model  are 
excellent  considering  that  experimental  variations  in  gas 
solubility  are  large.     (Errors  of  10.0%  are  expected.) 
Except  for  the  extremely  large  error  in  the  acetylene- 
acetone-n-hexane  system,   the  agreement  is  generally  good 
provided  the  excess  volume  is  included.     Otherwise,  the 
results  are  comparable  to  all  other  methods  which  are 
roughly  the  same.     Table  4-9  presents  average  deviations 
for  all  but  that  system.     To  indicate  the  difficulty  of 
describing  the  alcohol-water  mixtures   (22-40)  we  have 
separated  these  out.     It  is  clear  that  Mathias  model  pro- 
duces a  better  prediction  for  the  Henry  constant,  particu- 
larly in  aqueous  alcohols. 


117 


CN 

CN 

CN 

CN 

IN 

VD 

vo 

VD 

Oi 

O  1 

rH 

rH 

H 

r-l 

iH 

iH 

o 

00 

00 

CO 

00 

00 

CO 

c 

(N 

(N 

fNI 

CM 

CN 

(1) 

a\ 

a\ 

a\ 

<yi 

<Ti 

CTi 

C  >i 

c 

>1 

c 

>i 

d 

c 

C  >i 

u 

iH 

rH 

■H 

r-H 

(0 

X 

(0 

«5 

fC 

dJ 

tn 

m 

CO 

CO 

CO 

CO 

IW 

cn  > 

CO 

> 

CO 

> 

CO 

> 

CO 

> 

to  > 

x: 

x: 

xi 

x; 

Q)  0 

0) 

0 

(U 

0 

(U 

0 

0) 

0 

(1)  0 

4-1 

4-1 

4-) 

4J 

4-1 

4J  C 

4J 

c 

4J 

c 

4-1 

c 

4J 

c 

■P  c 

•H 

•H 

•H 

•H 

•H 

■H 

•H  (tS 

•H 

■H 

(0 

•H 

•H 

•H  (0 

E 

E 

g 

6 

e 

E 

Q)  > 

0) 

> 

0) 

> 

Q) 

> 

0) 

> 

Q)  > 

w 

to 

CO 

CO 

CO 

M 

H 

h:? 

H 

H 

fj  H 

O  (C 

4-1 

Z 


CO 

o 


E 

Q) 
4-1 
CO 

>1 

CO 


in 


in 


m 


n 


O 

u 


o 
u 


CN 


O 
U 


O 
U 


O 
U 


in 


O 

u 


CN 

u 


CN 

u 


00 


u 


CM 

u 


CM 

u 


in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

rH 

rH 

rH 

rH 

H 

H 

rH 

rH 

rH 

rH 

rH 

rH 

CO 

CO 

CO 

CO 

00 

00 

00 

n 

n 

cn 

n 

CTi 

cri 

CTi 

CTi 

rH 

CN 

rH 

CN 

CN 

CN 

CN 

CN 

CN 

CN 

CM 

CN 

CN 

CN 

CN 

CN 

o 

0) 

CO 

c 

c 

c 

0 

(CJ 

>1 

4-1 

X 

X 

Q) 

0) 

<u 

u 

(1) 

x; 

x; 

c 

< 

c 

0 

0 

•H 

1 

0) 

•H 

•H 

rH 

rH 

CP 

a) 

N 

u 

o 

0 

o 

c 

nj 

>1 

>1 

•H 

0) 

o 

c 

c 

0 

u 

c 

c 

H 

rH 

•H 

•H 

0 

•H 

•H 

c 

c 

m 

(0 

>1 

CO 

rH 

rH 

rH 

u 

4J 

4-1 

(0 

4J 

4J 

Si 

x: 

4-1 

-H 

•H 

4J 

Q) 

0) 

Oi 

CI. 

4J 

4J 

c 

CO 

c 

c 

•rl 

u 

u 

Q) 

0) 

0 

> 

•H 

-a 

< 

1 

< 
f 

1 

f 

< 
1 

m 

as 

K 

s 

s 

rH 

(U 

(1) 

0) 

o 

1 

Q) 

1 

0) 

1 

Q) 

1 

0) 

1 

0) 

1 

0) 

0 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

CO 

0 

o 

0) 

a) 

0) 

o 

0 

(1) 

Q) 

Q) 

N 

N 

d 

N 

4J 

4-1 

3 

3 

u 

H 

c 

c 

rH 

c 

(U 

Q) 

rH 

rH 

rH 

rH 

(0 

0 

CD 

0 

u 

o 

0 

0 

0 

0 

u 

Eh 

m 

m 

PQ 

<: 

< 

E-i 

EH 

Eh 

CN 


CN 


118 


00 

00 

00 

r- 

iH 

iH 

Ij 

Q) 

Q) 

w 

(0 

•0 

4h 

CO 

CD 

1 J 

vo 

VD 

0\ 
\V 

U/ 

1 

ON 

w 

r\ 
\J 

i-< 

i-J 

r^ 

rH 

rH 

,_) 

to 

tu 

lU 

Lj 

jj 

JJ 

m 

c 

B 

0) 

c 

<U 

0) 

0) 

0) 

>i 

•H 

•H 

•H 

•H 

•H 

IM 

(TJ  ^ 

iH 

rH 

rH 

(C 

(0 

(0 

rt  ^ 

(1) 

(0 

(0 

(d 

(0 

Q) 

OJ 

<u 

2 

Z 

z 

z  ^ 

Z 

•H 

+J 

4J 

+> 

+J 

-0 

-0 

T3 

1 

1 

1 

1 

u  a\ 

+J 

+j 

■P 

3 

3 

3 

3 

C 

C 

c 

C  <Ti 

C  CTl 

3  rH 

•H 

•H 

•H 

•H 

0 

0 

0 

0 

OJ 

0 

(1)  rH 

0)  H 

Pn  — 

!Z 

2 

z 

2 

CQ 

m 

m 

CQ 

CQ 

O 

o  nj 
Z  TJ 


CO 
IT) 

O 


e 

0) 
W 

CO 


in 


(Ti 


CTl 


Z 


"5r 


CM 

Z 


in 


(N 


IX) 


z 


rH 


o 

U 


00 


O 
U 


CM 

O 

u 


o 


CM 

o 
u 


rsi 


CN 
CM 


Sh 


in 

(N 


in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

in 

iH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

rH 

H 

rH 

rH 

rH 

00 

00 

00 

n 

n 

ro 

ro 

00 

00 

00 

00 

00 

CTl 

o> 

r- 

rH 

cri 

cri 

cn 

cri 

00 

CTl 

r- 

00 

CM 

CM 

CM 

CM 

rsi 

CM 

r>q 

CM 

CM 

CM 

CM 

CM 

CM 

0) 

u 

a 

0) 

■p 

rH 

nj 

■C 

(0 

rH 

0 

+1 

4-) 

4J 

0 

CU 

c 

u 

0) 

0) 

c 

c 

0 

u 

(U 

rH 

0 

Q) 

rH 

x; 

<u 

N 

0 

0 

m 

c 

4-> 

c 

c: 

c 

u 

Sh 

V4 

Sh 

H 

(U 

x; 

0) 

0) 

(U 

CU 

(U 

u 

0) 

cu 

1 

N 

4-> 

>1 

s 

N 

CQ 

x: 

4-1 

4-1 

4-) 

CU 

CU 

0 

0) 

G 

<D 

x: 

1 

c 

1 

4J 

(0 

03 

(0 

4J 

+J 

c 

(0 

c 

0) 

•H 

-p 

CU 

B 

CU 

s 

ns 

H 

(0 

oa 

Q 

w 

u 

CQ 

u 

S 

1 

1 

1 

> 

1 

X 

1 

1 

1 

0 

1 

0 

1 

rH 

rH 

rH 

1 

rH 

(1) 

Q) 

rH 

H 

rH 

m 

CU 

MH 

CU 

0 

0 

0 

rH 

H 

0 

c 

0 

0 

0 

0 

c 

0 

G 

c 

c 

c 

0 

0 

w 

0 

0 

c 

c 

c 

0 

0 

(d 

(0 

c 

CI 

-p 

rH 

0 

4-) 

0 

4J 

x; 

x: 

x: 

to 

0 

x; 

x: 

x: 

rH 

CU 

rH 

CU 

4J 

4J 

4-) 

x; 

XI 

o 

>1 

4-1 

HJ 

x: 

o 

XI 

o 

CU 

CU 

0) 

4J 

4J 

< 

o 

u 

w 

w 

u 

< 

u 

< 

S 

s 

s 

u 
< 


CM 


119 


U 

c 

u 

Q) 
IH 
CU 


0) 
(0 

PQ 

c 

E 

•H 

(C  ^ 
2  ^ 

C  a\ 
CD  >H 

m  — 


S-l 
CD 
(0 

m 
c 

(0 
£ 

■H 

rtJ  ^ 

2  ^ 

I  ^ 

0)  iH 
CQ  — 


(0 

c 

o 

Eh 


(0 

ro 
C 

d 
o 

EH 


r- 


OS 


CP 
(0 

c 
o 

Eh 


10 

c 
o 


(0 
(0 

c 
d 

o 

Eh 


CP 

to 

C 

0 
Eh 


(0 
(0 

c 
o 

Eh 


(0 
Cn 
(C 
C 

O 
Eh 


(0 

«3 
C 
3 

O 
Eh 


—  to 


to 

CP 
to 
C 

o 

Eh 


to 

to 

C 

O 


C 
lO 

a\ 
iC  ' 
CP-' 
to 
C  -H 
3  to 
^  ^ 
O  to 
Eh  fct; 


4-1 

o  to 

o  (0 
2 


in 


o 


Eh 


in 


00 
CM 


in 


cn 
o 
n 


in 


in 


in 


in 


in 


in 


in 


in 


in 


CTi 
<N 


fM 


in 


n 

CM 


in 


to 

4-) 

c 
> 

H 

o 


to 
I 

H 

o 
c 

to 

XI 

+J 
w 


0) 

-p 

to 

I 

0 

c 
to 
x: 
-p 

w 


0) 

to 
s 
I 

H 

o 
c 
to 

Xi 

+J 
s 


0) 

+J 
to 
s 

I 

cH 

0 

c 
to 
x: 

4-> 

s 


to 
I 

0 

c 
to 

Xi 

■p 

0) 


u 

0) 

■p 

to 

I 

o 
c 
to 
J5 
-p 
w 


^^ 
<u 
+j 
to 

I 

rH 
O 

c 

to 


u 

4-) 

to 

I 

o 

G 

to 

X! 
W 


•P 

to 

I 

iH 

o 
c 

10 

x: 
-p 

0) 

s 


0) 
4-> 

to 
s 

I 

iH 

o 
c 
to 
x; 
-p 

Q) 
S 


u 

<D 

to 

I 

o 
c 
to 
x; 
+j 
a 


u 

0) 

+) 
to 

I 

rH 
O 
C 

to 
x: 


-p 

10 

I 

H 
0 

c 
to 
x; 
-p 

s 


■p 

to 

I 

H 

o 
c 
to 
x: 
+) 

0) 

s 


to 


Sh 
< 


)H 
< 


CM        CN        CM  <N 

o    o    o  o 


CN  CN 

O  O 


CM 

2 


2 


CN 

2 


CM  E 

2  U 


U 


e 
-p 

(0 

>1 

CO 


CN 


00 
CN 


CTi 
CN 


o 

ro 


ro 


ro 


00 
ro 


ro 


in 
ro 


ro 


ro 


CO 

ro 


ro 


o 


120 


> 

n 

o 

CO 

on 

CN 

1 

\jj 

CN 

ID 

1 

lO 

o 

o 

IT) 

00 

o 

CO 

CN 

* 

1 

CO 

■ 

o 

o 

CO 

CN 

1 
1 

IN 

CO 

CN 

PO 

u  1 

CN 

1 

CO 

LO 

CN 

rH 

> 

CO 

U  1 

1 A 

u  ] 

CN 

* 

1 
1 

CN 

CN 

in 

ro 

o 

O 

CN 

rn 

^ 

CO 

VD 

CN 

CN 

iH 

o 

"<f 

CN 

00 

• 

• 

CN 

CN 

rH 

r^ 

CO 

H 

o 

ro 

rn 

tN 

• 

1 

O 

■H 

fH 

■H 

e 

0) 

o 

4J 

rH 

O 

CO 

(N 

>. 

CO 

UU 

0) 

rH 

c 

U3 

Hi 

lU 

1 

"H 

(0 

(1) 

H 

Si. 

XI 

+J 

rc 

Eh 

(0 

nj 

C 

•rH 

c 

cn 

•H 

6 

c 

(U 

4J 

0) 

CO 

E 

>i 

CO 

0) 

0 

+J 

(0 

•O 

0) 

c 

x; 

0 

■p 

ft 

CO 

0) 

> 

(0 

u 

Si 

0 

o 

a* 

<D 

C 

>^ 

0 

CO 

■1 

> 

0 

u 

CO 

CP 

c 

c 

•r1 

•H 

CHAPTER  5 

AN  APPLICATION  TO  LIQUID-LIQUID  EQUILIBRIA 


5 . 1  Introduction 

As  pointed  out  in  Chapter  4 ,  the  fluctuation  solution 
theory  is  successful  in  describing  the  gas-liquid  equi- 
libria.    A  more  severe  test  to  the  model  is  its  application 
to  liquid-liquid  equilibria.     For  the  latter,  no  model 
covers  all  aspects  of  two  or  more  fluid  phases  in  equi- 
librium.    This  chapter  describes  application  of  the  model 
of  Chapter  4  to  this  problem. 

In  a  binary  system,  there  are  two  pure  component 
choices  for  the  lower  limit  in  the  integral  equations  of 
the  fluctuation  solution  theory.     For  each  reference,  the 
model  contains  one  binary  parameter.     Due  to  the  imperfec- 
tions in  the  model,  the  calculated  results  and  parameters 
are  different  for  the  different  references.     The  agreement 
with  data,  while  good,  is  not  perfect.     In  order  to  improve 
the  situation  we  try  different  ways  to  use  the  equations 
with  two  references.     One  of  them,  the  multifluid  approach, 
combines  the  calculation  of  different  references.     The  ex- 
tension to  multicomponent  systems  is  straightforward. 

Section  5.2  presents  results  using  one  parameter  for 
each  binary  pair.     Section  5.3  describes  the  use  of  two 
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parameters  for  each  binary.     In  section  5.4,  we  present 
the  multifluid  approach  and  its  results. 

5 . 2    Calculations  of  Liquid-Liquid  Equilibria 
Using  One  Reference 

Two  different  classes  of  models  have  been  used  in 
correlating  liquid-liquid  equilibria: 

a)  Models  which  empirically  describe  the  distribution 
ratios.     Laddie  and  Degalasen  (1978)  present  a  good  s\ammary 
of  these  methods.     All  of  them  suffer  from  the  disadvantage 
of  working  for  only  ternary  systems  and  of  requiring  ex- 
perimental ternary  tie  lines  as  input  information. 

b)  Models  which  have  some  semi theoretical  basis. 
S0rensen  et  al.    (1977)  present  an  excellent  review  with 
extensive  comparisons  between  methods.     In  general,  these 
models  describe  the  excess  Gibbs  free  energy  as  a  function 
of  composition  and  temperature.     The  UNIQUAC   (Abrams  and 
Prausnitz,  1975)  equation  is  the  most  convenient  for  cor- 
relating data,  but  it  is  not  completely  successful. 

Most  binary  liquid  mixtures  have  an  upper  critical 
solution  temperature    due   to    the   balance  of  the  entropy 
of  mixing  against  the  interaction  energy  of  the  molecules. 
More  complex  mixtures,  such  as  amines  in  solution  with 
water,  generate  temperature-composition  phase  diagrams 
with  an  upper  and  a  lower  critical  solution  temperature. 
This  behavior  can  be  understood  by  the  presence  of  highly 
directional  interactions ,   such  as  hydrogen  bonds .  The 
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models  mentioned  above  do  not  provide  a  good  correlation 
of  these  closed-loop  diagrams.     For  these  cases,  the  most 
successful  theories  are  applications  of  the  Ising  model 
(gas-lattice  models   (VJheeler  and  Andersen,  1980;  Goldstein 
and  Walker,  1983) ) . 

In  order  to  predict  binary  and  multicomponent  liquid- 
liquid  equilibria,  it  is  necessary  to  have  a  method  for 
predicting  activity  coefficients.     The  UNIFAC  (Fredenslund 
et  al.,  1975)  method,  a  group  contribution  model,  is  one 
of  such  methods.     Magnussen  et  al.    (1981)  built  a  group 
interaction  parameter  table  suited  for  liquid-liquid 
equilibria.     Their  results  indicate  that  only  a  semi- 
quantitative prediction  is  possible. 

We  shall  use  the  model  for  the  direct  correlation 
function  integrals  described  in  section  3.2,  together  with 
equations  given  in  section  2.4,  to  correlate  liquid-liquid 
equilibria  in  binary  and  multicomponent  mixtures. 

The  procedure  to  do  the  calculations  involves  an 
iterative  calculation  with  the  following  six  steps: 

1)  Fitting  of  pure  component  characteristic  parameters. 
We  evaluate  T*  and  V*  with  equation  (4-4)  and  P-V-T  data. 
Table  5-1  presents  the  characteristic  parameters  for 
substances  involved  in  our  liquid-liquid  equilibrium  calcu- 
lations . 

2)  Calculation  of  final  solution  density.  A  numerical 
solution  is  found  for  equation  (4-12) .  To  do  this,  we  need 
guesses  of  the  binary  parameter  and  mixture  composition. 
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3)  Calculation  of  activity  coefficients.  Activity 
coefficients  are  found  from  equation  (A-4)  using  the 
binary  parameter,  composition,  and  calculated  solution 
density  of  step  2 . 

4)  Confirmation  of  equilibrium  composition.     Steps  2 
and  3  are  done  for  each  phase  present  in  the  equilibrium. 
The  isoactivity  condition  expressed  by  equation  (2-50)  is 
now  checked.     If  the  calculated  activity  coefficients  do 
not  satisfy  equation  (2-50) ,  a  new  composition  of  the 
solution  in  each  phase  is  guessed,  and  steps  2  and  4  are 
repeated . 

5)  Confirmation  of  stability.     If  the  two  phases 
follow  the  isoactivity  requirement,  we  check  material 
stability.     If  the  phases  do  not  satisfy  this  condition, 
a  new  composition  is  guessed,  and  steps  2  and  5  are  re- 
peated . 

6)  Comparison  of  experimental  and  calculated  concen- 
trations.    If  the  error  is  larger  than  expected,  the  binary 
parameters  are  changed,  and  steps  2  to  6  are  repeated 
until  further  changes  do  not  improve  the  agreement. 

Winnick  and  Powers  ( 1966a, b)  made  a  detailed  study 
of  the  acetone-carbon  disulphide  system.     At  0°C  and 
atmospheric  pressure  this  system  is  completely  miscible, 
but  at  pressures  over  5000  atmospheres  partial  miscibility 
can  be  observed.     They  measured  the  mixing  vol;ame  as  a 
function  of  pressure  and  temperature,  and  the  phase  be- 
havior at  high  pressure.     They  used  the  volumetric  data 
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to  calculate  equilibrium  composition.     Figure  5-1  shows 
the  observed  phase  behavior  that  calculated  with  the  volu- 
metric data  and  that  calculated  using  the  fluctuation 
solution  theory.     Agreement  between  calculated  and  experi- 
mental is  better  using  direct  correlation  function  integrals 
than  using  volumetric  data.     Attempts  to  improve  the  agree- 
ment by  changing  the  values  of  the  characteristic  parameters 

(T*  and  V*)  were  unsuccessful  and  Figure  5-1  gives  the 
best  representation  of  the  method.     However,  the  uncer- 
tainties in  this  system  are  large  since  even  the  calcu- 
lated phase  behavior  of  Winnick  and  Powers  does  not  agree 
with  experiment. 

An  alternative  system  is  that  of  Keyes  and  Hildebrand 

(1917)  who  studied  the  aniline-hexane  system.     They  measured 
solubility,  mixing  volume,  heat  of  mixing,    surface  tension, 
and  compressibility  in  the  mixture.     Figure  5-2  displays 
the  experimental  points  and  the  calculated  phase  envelope 
from  the  present  method.     The  agreement  is  quite  good.  To 
observe  the  effect  of  the  characteristic  parameters  T*  and 
V*  on  the  phase  envelope,  we  changed  their  values  and  made 
calculations  using  the  six  step  procedure.     These  changes 
do  not  affect  the  shape  of  the  phase  envelope.  Rather 
they  move  its  position  in  a  temperature-composition  diagram. 
Figure  5-3  presents  a  schematic  picture  of  these  displace- 
ments using  the  upper  critical  solution  temperature  (CST) . 
These  changes  do  not  produce  an  improved  correlation  of 
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Figure  5-1.     Pressure  dependence  of  liquid-liquid 
equilibria  for  acetone   (1) -carbon 
disulphide   (2)   at  T  =  273.15  K. 
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Figure  5-2.     Temperature  dependence  of  liquid- 
■  liquid  equilibria  for  aniline  (1)- 
hexane   (2)   at  1  atm. 
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Movement  of  upper  critical  solution 
temperature  and  liquid-liquid  phase 
envelope  of  aniline   (l)-n-hexane  (2) 
system  with  changes  in  parameters. 
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the  liquid-liquid  equilibria  for  aniline-hexane.  The  phase 
envelope  in  Figure  5-2  is  our  best  calculation. 

In  order  to  compare  the  results  from  the  fluctuation 
solution  theory  with  those  of  UNIQUAC,  an  analysis  of  the 
aniline-cyclohexane  is  made.     Angelescu  and  Guisca  (1942) 
measured  the  solubility  in  aniline-cyclohexane  mixtures. 
Their  results  are  presented  in  Table  5-2.     Eckert  et  al. 
(1980)   examined  this  system  with  the  UNIQUAC  model.  Table 
5-2  presents  their  findings.     They  analyzed  an  extension 
of  UNIQUAC  using  four  parameters  by  allowing  the  original 
parameters  to  be  linear  functions  of  temperature.  Though 
this  provides  a  better  correlation  than  the  original  two- 
parameter  model,  we  have  chosen  to  do  the  comparison  with 
the  original  UNIQUAC.     Table  5-2  presents  the  composition 
calculated  with  the  fluctuation  solution  theory  using  pure 
aniline  as  the  reference  state.     Figure  5-4  displays 
comparisons  between  experimental  and  calculated  composi- 
tions.    For  this  particular  system,  the  correlation  given 
by  the  fluctuation  solution  theory  is  better  than  the  one 
given  by  the  original  UNIQUAC  model. 

5.3     Calculations  of  Liquid-Liquid  Equilibria 
Using  Two  References 

The  main  problem  with  the  model  described  in  section 
5.2  is  that  it  contains  only  one  parameter  whose  value 
depends  on  the  reference  state.     Thus  the  compositions  in 
the  two  phases  in  a  binary  system  cannot  be  matched.  An- 
other way  to  look  at  the  problem  of  the  model  is  by  doing 
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Figure  5-4.     Temperature  dependence  of  liquid-liquid 
equilibria  for  aniline   ( 1) -cyclohexane 
(2)   at  P  =  1  atm. 
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density  calculations  using  different  references.  Figure 
5-5  displays  the  differences  that  exist  when  a  different 
reference  is  chosen  to  do  the  calculation  of  the  solution 
density.     The  density  of  one  pure  component  cannot  be 
predicted  accurately  using  only  the  reference  density  and 
equation   (4-12) . 

Experience  shows  that  in  predicting  ternary  systems, 
the  best  starting  point  is  a  perfect  match  for  the  immisci 
ble  pairs.     Other  models  have  overcome  this  difficulty  by 
having  two  binary  parameters.     Thus,  we  shall  try  two 
alternatives  to  introduce  a  second  binary  parameter  in 
our  model  and  provide  similar  matching: 

a)  Use  one  binary  parameter  for  each  phase  with  a 
single  reference  for  both  phases. 

b)  Use  of  different  references  for  each  phase,  each 
having  its  own  binary  parameter.     For  the  case  of  Figure 
5-5  this  alternative  implies  using  hexane  as  reference  in 
the  hexane-rich  side  and  aniline  in  the  aniline-rich  side. 

Table  5-3  presents  the  calculations  for  the  system 
aniline-cyclohexane  at  10°C  using  alternatives  a  and  b. 
Two  observations  can  be  made: 

1)  One  reference,  two  parameters  (alternative  a)  is 
only  as  good  as  one  reference  one  parameter. 

2)  Two  references,  two  parameters   (alternative  b)  is 
equivalent  to  one  reference,  one  parameter  in  each  phase. 
Most  importantly,  calculations  of  the  binodal  curve 
yielded  poor  or  no  convergence. 
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Figure  5-5.     Solution  volumes  for  aniline  (1)- 
hexane   (2)   at  T  =  298.15  K. 
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Table  5-3 

Comparison  of  Tie  Line  Calculations  Using  Different 
Binary  Parameters  and  Different  References  for 
Aniline  (1) -Cyclohexane  (2)  at  10°C 


Reference:  Cyclohexane 

a.  Binary  parameter  =  .026 

x^,  phase'  =  .1118  x^ ,  phase"  =  .9121 

b.  Binary  parameter  =  .03 

x^,  phase  '  =  .0666  -x.^,  phase"  =  .9493 

c.  Binary  parameter  cyclohexane  side  =  .026, 
binary  parameter  aniline  side  =  .03 

x^,  phase'  =  .0621  y^^t  phase"  =  .9491 

d.  Binary  parameter  cyclohexane  side  =  .03, 
binary  parameter    aniline   side  =  .26 

x^,  phase'  =  .1402  X2,  phase"  =  .9133 

Reference:  Aniline 

a.  Binary  parameter  =  .026 

x^,  phase'  =  .0786  x^,  phase"  =  .9401 

b.  Binary  parameter  =  .03 

x^,  phase'  =  .0464  x^^,  phase"  =  .9657 

c.  Binary  parameter  cyclohexane  side  =  .026, 
binary  parameter  aniline  side  =  .03 

x^,  phase'  =  .0783  x^,  phase"  =  .9333 

d.  Binary  parameter  cyclohexane  side  =  .03, 
binary  parameter  aniline  side  =  .026 

x^,  phase'  =  .0464  x^,  phase"  =  .9672 

Reference:     Aniline  and  cyclohexane 

a.  Binary  parameter  cyclohexane  side  =  .026, 
binary  parameter    aniline   side  =  .03 

Xj^,  phase'  =  .1215  x^^,  phase"  =  .9679 

b.  Binary  parameter  cyclohexane  side  =  .03, 
binary  parameter  aniline  side  =  .026 

x^,  phase'  =  .0659  x^,  phase"  =  .9392 

Experimental   (Angelescu  and  Giusca,  1942) . 

Xj^,  phase'  =  .  072  x^^ ,  phase"  =  .823 


136 

In  order  to  clarify  the  results  obtained  in  the  calcu- 
lations using  two  binary  parameters  and  two  references,  we 
look  at  the  shape  of  the  chemical  potential  against  com- 
position.    The  system  tetraf luormethane-methane  is  appro- 
priate because  it  shows  partial  miscibility  between  87  K 
and  95  K  at  atmospheric  pressures,  and  there  is  extensive 
information  about  the  system  (Croll  and  Scott,  1958;  Thorp 
and  Scott,  1956;  Simon  and  Knobler,  1971;  Paas  and  Schneider, 
1979) .     The  results  of  the  chemical  potential  for  the 
tetraf luormethane-methane  mixtures  are  presented  in  Tables 
5-4,  5-5,  and  5-6.     In  Tables  5-4  and  5-6  we  can  see  that 
an  increase  in  the  binary  parameter  produces  an  increase 
in  the  chemical  potential.     The  increase  is  larger  at 
infinite  dilution  for  both  components;  near  pure  composi- 
tion, the  variation  is  small.     Table  5-5  shows  the  obvious 
result  that  the  chemical  potentials  are  equal  when  the 
references  and  the  binary  parameters  are  equal.     While  the 
chemical  potential  curve  for  a  particular  component  is 
completely  different  when  calculated  with  two  references 
and  two  parameters  near  pure  component,  it  is  practically 
the  same  as  when  calculated  with  one  reference  and  two 
parameters . 

Now,  it  is  possible  to  explain  the  behavior  of  the 
calculations  in  the  aniline-cyclohexane  system.     We  made 
the  observation  that  the  alternative  a  (use  of  one  reference 
and  two  paramters)   is  equivalent  to  one  reference,  one 
parameter.     The  equilibrium  is  found  by  matching  two 
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chemical  potential  curves.     In  these  cases,  the  curves 
are  identical  at  high  concentrations  of  the  component. 
As  a  consequence,  they  cannot  be  distinguished.     The  be- 
havior of  alternative  b  is  due  to  the  same  effect;  each 
phase  follows  the  behavior  of  the  curves  near  pure  com- 
ponent for  chemical  potential  matching. 

5 . 4     Calculations  of  Liquid-Liquid  Equilibria 
Using  a  Multifluid  Approach 

To  overcome  the  limitation  of  using  one  parameter,  we 
explore  a  multifluid  approach  that  gives  for  a  binary 
system 


where  the  supraindex  on  the  left  indicates  the  substance 
chosen  as  reference  to  do  the  calculations   (e.g.,  ^Yj^  is 
the  activity  coefficient  of  component  i  calculated  using 
component  1  as  reference)  .     By  definition  of  In        we  can 
write 


In  Yi 


^1  "^i      ^2  "^i 


(5-1) 


n 


RT    Z    X.  In 
1  1 


(5-2) 


then,  combining  equations   (5-1)   and  (5-2) 


+  X 


(5-3) 
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Excess  enthalpy  and  excess  voliune  are  related  with  excess 
Gibbs  energy  through 


8G 


E 


T,x 


(5-4) 


3T  P,x 


SI 
■  T 


(5-5) 


Now,  using  equation   (5-3)   in  equations   (5-4)   and  (5-5) 


=    x^  V  +  X2  V 


(5-6) 


and 


=     x^  +  X2 


(5-7) 


Table  5-7  shows  the  results  of  the  tetraf luormethane- 

methane  system  in  the  partially  and  completely  miscible 

1  2 

regions.     The  binary  parameters  '^ll  obtained 

from  liquid-liquid  equilibria  at  89  K  and  estimated  from 
the  excess  Gibbs  free  energy  at  110  K.     A  linear  function 
of  temperature  was  used  for  all  the  other  data.     The  agree- 
ment between  calculated  excess  Gibbs  energy  and  excess 
volume  with  experimental  data  is  quite  good. 

The  aniline-cyclohexane  system  in  the  partially  and 
completely  miscible  regions  is  examined  thoroughly  here. 
The  calculations  cover  liquid-liquid  equilibrium  composition 
excess  Gibbs  free  energy,  excess  volume,  and  excess  enthalpy. 
In  Table  5-8  the  binary  parameters    k, ~  and    k, ^  were 
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obtained  from  liquid-liquid  equilibria  at  278.15  K, 
288.15  K,  and  298.15  K.     A  linear  interpolation  was  used 
for  all  the  other  calculations.     In  Table  5-9  the  linear 
interpolation  for  the  binary  constants  was  done  with  the 
liquid-liquid  equilibria  at  278.16  K  and  298.16  K.  Figure 
5-6  displays  the  results  of  Table  5-9.     In  Table  5-10  the 
binary  parameter  function  of  temperature  was  interpolated 
between  the  278.15  K  and  392.46  K  data.     Table  5-11  presents 
the  results  of  calculations  with  UNIQUAC  model  using 
parameters  given  in  the  literature  (Nicolaides  and  Eckert, 
1978a) .     Finally,  Table  5-12  compares  experimental  and 
calculated  results  using  the  UNIQUAC  model  and  the  fluctua- 
tion solution  theory. 

From  all  these  results  we  can  infer  the  following  con- 
clusions for  fluctuation  solution  theory: 

a)   In  the  correlation  mode  with  a  linear  dependence  in 
the  binary  parameters,  the  description  of  the  liquid-liquid 
equilibria  is  reasonably  good  (Table  5-8)  with  good  results 
also  obtained  in  the  completely  miscible  region.  Results 
of  the  same  quality  are  obtained  with  only  two  sets  of 
liquid-liquid  equilibrium  data  (Table  5-9)  or  with  one  set 
of  data  in  the  partially  miscible  region  and  another  in 
the  miscible  region  (Table  5-10)  .     The  results  calculated 
from  the  fluctuation  solution  theory  are  better  than  the 
results  calculated  using  UNIQUAC  (Table  5-11)  which  has 
large  errors  in  the  miscible  region. 
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Figure  5-6.     Correlation  of  liquid-liquid  equilibria 
and  excess  Gibbs  energy  for  aniline  (1) 
and  cyclohexane   (2)   using  multifluid 
references . 
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b)  The  value  of  the  excess  enthalpy  is  of  the  order 
of  experiment,  but  its  dependence  on  temperature  is  wrong 
(Table  5-12) . 

The  calculation  of  properties  for  the  water-butanol 
system  is  a  more  severe  test  for  any  model,  given  the 
particular  characteristic  of  water  (hydrogen  bonding, 
polarity,  etc.).     Table  5-13  shows  the  results  for  the 
water-butanol  system.     The  agreement  between  experimental 
and  calculated  composition  is  not  as  good  as  aniline- 
cyclohexane,  but  it  is  not  too  bad.     In  the  water-butanol 
system  the  binary  parameters  do  not  have  a  linear  varia- 
tion with  temperature   (Figure  5-7) .     Figure  5-8  shows  the 
water  activity  at  353.16  K  and  at  393.16  K.     The  flatness 
of  the  curves  at  high  water  concentration  explains  the 
anomalous  results  of  Table  5-13;  a  very  small  change  in 
water  activity  produces  a  large  change  in  the  water  mole 
fraction . 

The  extension  of  equation  (5-1)   to  an  n-component 
system  is 

n 

In  Yi  =     I     X.   In  \.  (5-8) 
^       j=l     3  1 

In  equation  (5-8)   the  -^y^^  depends  on  binary  parameters 
through  equation  (A-  4 ) .     The  multicomponent  system  brings 
a  "^k^j  parameter  where  l^^i^j  .     From  binary  data  we  calcu- 
late ■'"kj^j  and  ^k^ . ,  and  "^k^ .  is  obtained  using 
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Table  5-13 

Comparison  of  Liquid-Liquid  Compositions  for 
Water  (l)-Butanol  (2) 


T,  K 

^1' 

phase ' 

^1' 

phase 

Exp. 

Calc. 

Exp. 

Calc . 

313.15 

.9831 

.9847 

.528 

.514 

323.15 

.9835 

.9839 

.543 

.538 

333.15 

.9833 

.9833 

.560 

.560 

343.15 

.9827 

.9827 

.581 

.581 

353.15 

.9823 

.9402 

.597 

.606 

363.15 

.9796 

.9397 

.639 

.623 

373.15 

.9760 

.9391 

.677 

.640 

383.15 

.9720 

.9806 

.712 

.646 

388.15 

.9660 

.741 

393.15 

.9550 

.9379 

.778 

.671 

The  experimental  data  are  from  Hill  (1928)  . 
The  calculated  data  are  from  fluctuation  solution 
theory . 
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— '  1  1  I  I  L 

293  313  333  353  37l  TO 

Temperature,  K 


Figure  5-7.     Binary  parameters  in  equation   (5-1)  for 
water   (1)   and  butanol   (2)   at  P  =  1  atm. 


153 


1.06 

X 

K 

X 

X 

X 

X 

X 

1.04 

X 

X 

—      X  "w^^ 

.98 

T  =  353.16  K 

X 

X  _ 

Q  A 

.• 

K 

T  =  393.16  K 

X 

.9 

! — «  1 — 1  1 

XI         1         1  , 

•5     .6     .7     .8     .9     1.0  .5     .6     .7     .8     .9  1.0 


Figure  5-8.     Water   (1)   activity  in  butanol  at 
P  =  1  atm. 
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\.  =  i  (\.   +  \.)  (5-9) 

Table  5-14  and  Figure  5-9  show  the  results  of  equi- 
librium calculations    in  n-heptane-cyclohexane-aniline 
system.     The  numbers  at  low  concentration  of  cyclohexane 
are  quite  good,  but,  as  noted  in  Table  5-14,  in  the  mole 
heptane  fraction  range  from  0.6  to  3.5  no  solution  was 
found.     This  occurs  because,  as  Appendix  C  shows,  the  multi- 
fluid  approach  of  equation  (5-8)  does  not  satisfy  the 
Gibbs-Duhem  equation. 

Returning  to  the  one-parameter,  one-reference  model. 
Table  5-15  and  Figure  5-10  present  the  results    of  calcu- 
lations on  n-heptane-aniline-cyclohexane  except  that  the 
pure  component  parameters  T*  and  V*  were  adjusted  so  that 
the  calculation  is  independent  of  the  reference.  The 
characteristic  temperature  of  n-heptane  was  changed  to  fit 
the  experimental  n-heptane  volume  for  the  aniline  ref- 
erence.    Also  the  characteristic  temperature  of  cyclohexane 
was  adjusted  to  fit  the  pure  cyclohexane  volume  with  the 
aniline  reference  .    The  characteristic  temperature  values 
are  shown  in  Table  5-15.     The  numbers  obtained  with  this 
new  type  of  approach  are  not  good,  particularly  on  the 
aniline  side,  as  shown  in  Figure  5-10  and  Table  5-15. 

The  limited  testing  of  the  liquid-liquid  calculations 
method  using  fluctuation  solution  theory  shows  some  weak- 
nesses   in   the    perturbed  hard-sphere  models  presented. 
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Cyclohexane 


Figure  5-9.     Liquid-liquid  equilibria  in  the  system 
heptane   ( 1) -cyclohexane   (2) -aniline  (3) 
at  T  =  298.15  K  using  multifluid  references 
and  two  parameters  per  binary 
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Cyclohexane 


Aniline  0.2  0.4  0.6  0.8  Heptane 


Figure  5-10.     Liquid-liquid  equilibria  in  the  system 
heptane   (1) -cyclohexane   (2) -aniline  (3) 
at  T  =  298.15  K  using  single  fluid  ref- 
erences and  one  parameter  per  binary. 
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In  general  for  binary  systems,  the  results  are  ac- 
ceptable.    Not  only  does  the  model,  using  the  multifluid 
approach,  correlate  the  phase  envelope  well,  it  provides 
multiproperty  correlation  with  linear  temperature  parameters. 
On  the  other  hand,  the  UNIQUAC  model  does  not  show  these 
features.     This  comparison  between  models  indicates  a  better 
temperature  dependence  built  in  the  perturbed  hard  sphere 
model.     We  mentioned  only  the  existence  of  closed-loop 
diagrams  in  binary  systems.     These  mixtures  have  an  upper 
and  a  lower  critical  solution  temperature  due  to  the  bal- 
ance of  two  opposite  effects  that  have  different  temperature 
dependence.     It  is  clear  that  the  models  presented  in  this 
chapter  with  linear  dependent  parameters  cannot  correlate 
the  closed-loop  system.     We  have  made  calculations  using 
the  multifluid  approach  that  confirm  this  limitation.  As 
a  result,  we  cannot  present  modifications  to  the  model  to 
improve  its  temperature  dependence. 

In  ternary  systems,  the  results  are  not  good.  The 
multifluid  approach  violates  the  Gibbs-Duhem  equation,  and 
the  mixing  Gibbs  energy  surface  generated  by  these  activity 
coefficients  does  not  possess,  at  some  compositions,  a 
common  tangent  plane  in  the  regions  corresponding  to  two 
different  phases.     The  use  of  a  one-parameter  model  is  not 
good  even  for  ternary  systems  with  two  inmiscible  binaries 
which  should  be  the  easiest  to  correlate. 

The  results  in  ternary  systems  indicate  that  the 
composition  dependence  in  the  model  is  not  quite  right. 
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and  this  imperfection  is  a  very  important  weakness.  In 
the  next  chapter  we  make  some  suggestions  to  improve  the 
composition  dependence  of  direct  correlation  function 
integral  models. 


CHAPTER  6 

MOLECULAR  THERMODYNAMIC  MODELS  FOR  DIRECT  CORRELATION 

FUNCTION  INTEGRALS 

6 . 1  Introduction 

In  applying  fluctuation  solution  theory  to  liquid- 
liquid  equilibria,  the  methods  used  up  to  this  point  have 
proven  inadequate.     The  one-parameter  model  for  binary  sys- 
tems in  Section  5.2  is  reasonably  good,  but  it  is  not  ade- 
quate for  ternary  mixtures .    Alternatives  for  overcoming 
the  one-parameter  difficulty  also  had  problems  that  led 
to  their  rejection.     Further,  in  Section  4.3  it    is  indi- 
cated that  the  hard  sphere  mixture  may  not  be  a  good 
reference  for  describing  direct  correlation  function 
integrals . 

An  additional  problem  appears  in  the  infinite  dilution 
behavior  of  the  direct  correlation  function  integrals .  Even 
in  binary  systems  data  show  that  the  C^^  depend  on  unlike 
interactions  between  the  components.     The  basic  model,  or 
its  modifications,    do   not    take    this  into  account,  and 
the  only  interaction  that  appears  at  infinite  dilution  is 
that  of  hard  spheres .     One  way  to  get  rid  of  this  limita- 
tion is  to  include  another  term  in  the  expansion,   such  as 
quadratic  in  composition  and  density.     Another  way  is  to 
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look  at  different  expressions  that  are  obtained  from  models 
for  classical  excess  Gibbs  energy  expressions. 

Section  6.2  analyzes  hard  sphere  mixtures  as  the  ref- 
erence.    Section  6.3  presents  the  idea  of  including  an 
extra  term  in  the  expansion .     Section  6 . 4  examines  the 
excess  Gibbs  energy  models. 

6 . 2    Behavior  of  Perturbed  Hard  Sphere 

In  Chapter  4  we  concluded  that  for  mixtures  the 
particular  choice  of  hard  sphere  equation  for  the  reference 
term  in  equation  (4-1)   is  not  crucial.    We  now  address  the 
question  of  whether  the  hard  sphere  expression  plus  a 
linear  density  correction  term  is  really  a  good  way  to 
model  direct  correlation  function  integrals. 

In  this  section,  we  work  with  equations  very  similar 
to  those  of  Chapter  4,  except  with  different  parameteri- 
zations .     The  characteristic  temperature  and  characteristic 
volume  were  the  basic  parameters  in  Section  4.2.  Here, 
they  are  replaced  by  the  hard  sphere  diameters  of  the  pure 
components.     The  direct  correlation  function  integrals  at 
constant  temperature  can  be  written 

C^j(x,p)  =  C^j^^(x,p,a^,a2)       P  ®ij  ^^"^^ 

where  x  is  the  mole  fractions  vector,  p  is  the  density  of 
the  solution,   a.   is  the  hard  sphere  diameter  of  component 
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i,  and  B. .  is  the  coefficient  of  the  correction  term.  In 
ID 

addition 


^12  "  ^12   ^^"^^  ^^"^^ 


with  3  as  a  binary  parameter  used  to  match  the  calculated 
and  the  experimental  activity  coefficients . 

The  values  of  the  hard  sphere  diameters  are  changed 
at  constant  density  to  make  the  B  term  a  constant  while 
simultaneously  matching  the  experimental  direct  correlation 
function  integrals . 

Table  6-1  presents  the  results  for  the  acetone  (1)- 
carbon  disulphide  (2)   system.    The  column  under  shows 
the  values  obtained  using  equation  (6-1)  with  B^^  taken  as 
the  average  of  the  difference  between  the  experimental  C^^ 
and  the  hard  sphere  C^j    (B^^).     At  constant  temperature  the 
coefficients  of  the  correction  term  and  the  hard  sphere 
diameters  are  constant,  so  changing  the  density  of  the 
mixture  should  not  change  their  values . 

Table  6-2  presents  the  same  kind  of  results  as  Table 
6-1  but  at  a  different  mixture  density.     Unfortunately  the 
average  B^^  has  increased  around  10.0%  rather  than  being  the 
same  as  in  Table  6-1. 

Table  6-3  shows  a  comparison  between  C^^  from  differ- 
ent sources:     experimental,  calculated  using  equation  (4-6) 
with  characteristic  parameters  from  Table  5-1,  and  k^2  = 
.04,  calculated  using  equation  (6-1)  with  B  =  .0  and 
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calculated  using  equation  (6-1)  with  6  =  -.035.     The  value 
of  g  =  -.035  was  obtained  matching  the  experimental  and 
the  calculated  activity  coefficients.     The  number  in 
parentheses  in  each  column  are  the  difference  with  experi- 
ment. 

As  pointed  out  before,   the  difference  with  experiment 
increases  when  C^^  is  calculated  (using  equation  (6-1))  at 
a  density  different  from  that  in  which  B^^  was  evaluated. 
We  can  make  a  similar  observation  for  the  C^^^  calculated 
with  equation  (4-6) ;  the  difference  with  experiment  increases 
with  increasing  density.     Another  observation  is  that  the 
*^11  ^^^'^^'^        match  the  activity  coefficients  are  very 
different  from  the  experimental  C^^^. 

Figure  6-1  displays  the  mixture  volume  against  mole 
fraction.     The  figure  compares  the  experimental  curve  with 
results  obtained  using  direct  correlation  function  inte- 
grals calculated  with  equation   (4-6)   and  calculated  with 
equation   (6-1) .     While  the  results  using  different  refer- 
ences are  somewhat  different,  the  major  feature  is  the 
obvious  disagreement  with  experiment. 

Table  6-4  presents  results  of  the  volumes  and  the 
activity  coefficients  using  direct  correlation  function 
integrals  calculated  with  equation   (6-1)   using  B^^  as  the 
average  value    of   Table  6-1.     Results  with  equation  (4-6) 
for  the  are  also  included  in  the  table.     Small  changes 

in  the  value  of  the  coefficients  of  the  linear  density  term 
produce  large  changes  in  the  activity  coefficient  but 
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 J  ,  1 

Data  of  Winnick  and  Powers  (1969a) 

Characteristic  parameters  from  pure 
compressibilities 

Calculated  using  equation  (6-1)  with 


I  I  I  I  I 

0.0  0.25  0.50  0.75         x,  1.0 


Figure  6-1.     Mixture  volume  calculations  for  acetone 
(1) -carbon  disulphide   (2)   at  273.15  K 
and  1  atm. 
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only  a  small  change  in  the  mixture  volume.  Changing 
through  the  binary  parameter  produces  no  change  in  the  cal- 
culation of  the  pure  component  volume  and  only  a  slight 
change  at  other  concentrations.     A  change  in  B^^  produces 
a  small  change  at  all  concentrations .     The  changes  in  ac- 
tivity coefficient  with  a  change  in  either  B^^  or  B^^ 
very  large.     Table  6-5  compares  the  difference  in  the  excess 
free  Gibbs  energy  and  in  the  activity  coefficient  of  acetone 
calculated  using  acetone  as  the  reference.     The  changes 
here  are  more  significant;  a  change  of  1.0%  in  B^^^  produces 
a  change  of  20.0%  in  the  excess  free  Gibbs  energy  for  an 
equimolar  mixture  and  a  change  of  more  than  100.0%  in 
the  activity  coefficient  at  infinite  dilution.     The  changes 
produced  by  adjusting  B^2  to  match  the  activity  coefficient 
are  even  larger . 

We  can  now  say  that  the  hard  sphere  reference  is  not 
generally  a  good  concept.     First,  the  coefficients  of  the 

correction  term  must  be  density  dependent  to  correlate  the 

hs 

data.     Second,  because  of  the  difference  in  sign  of  C. . 
and  the  correction  term,   the  calculated  thermodynamic 
properties  are  very  sensitive  to  values  of  these  coeffi- 
cients . 


6 . 3    A  Three  Term  Expansion  for  the  Direct 
Correlation  Function  Integrals 


In  Chapter  3  we  showed  the  experimental  behavior  of 
nonideal  mixtures.     We  also  found  an  expression  for  the 
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direct  correlation  function  integrals  as  functions  of  the 
mixture  compressibility,  partial  molar  volume,  and  deriva- 
tives of  the  activity  coefficients.     The  particular  case  of 
C^^  for  a  binary  1-2  can  be  expressed  as 

1 


col        ^^1  ^                d  m  Yi 
(1-C-.)    ^  =  — ^-   +    —5  i 

^1  V2°<2°  RT  ^^1 


(6-3) 


where  the  superscript  <»1  indicates  that  the  mole  fraction 
of  component  one  goes  to  zero.     From  equation  (6-1)  we  can 
get  the  infinite  dilution  expression  for  a  perturbed  hard 
sphere  equation  for  C^^ 

(1-C^^)"^  =   [(1-C^^)^^]"^  +  B^3_  (6-4) 

Comparing  equation  (6-3)  with  equation   (6-4) ,  we  can  point 
out  an  important  weakness  in  the  perturbation  model.  Equa- 
tion (6-3)    is  highly  dependent  on  1-2  interaction,  but  in 
equation  (6-4)   the  interaction  1-2  is  only  in  the  hard 
sphere  mixture  term.     This  makes  explicit  the  inadequacy 
of  the  linear  perturbation  term. 

Additional  terms  can  be  added  to  equation  (6-1)  to 
take  account  of  this  defect.     These  additional  terms  must 
have  thermodynamic  consistency 
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where  P  stands  for  the  pressure  of  the  system.  As 


=x,  (1-C,J   +  x^(l-C,^)  (6-6) 


3p,         1'"  "11'       "2^"  "12' 


||-  =  x^(l-C^2)   4-  X2(l-C22)  (6-7) 


Replacing  equation  (6-6)   and  (6-7)   in  (6-5)  we  get 


^1   fP^T  (1-^11^   -  P  ih  ^"-^12^   ^  ^11  -  ^12^ 


=  ^2   tP  (1-^22)   -  P  (1-^12)   ^  S2  -  ^12^ 


(6-8) 

The  only  kind  of  extra  term  that  satisfies  equation  (6-8) 
is  of  the  third  virial  type.  The  expression  for  one  such 
alternative  is 


C^.    (x,£)   =       .  ^(x,p,  0^,02)  +pB. 


The  infinite  dilution  limit  of  equation  (6-9)  is  for 
i  =  j  =  1 


1  ,001  2 

(l-C^^)"-^  =   [(1-C^^)^^]       +  B^^  +(P2°)  D^^2  (6-10) 
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where  the  term  D 


112 


contains  information  about  the  1-2 


interaction.     Reordering  equation  (6-9)  gives 


hs 


2 


k  ijk 


hs 


+  p  b:  . 


(6-11) 


where  Bl.   is  the  coefficient  of  the  correction  term  that  is 


linearly  density  dependent. 

Unfortunately,  the  results  of  Table  6-3  and  the  plots 
in  Chapter  4   (e.g..  Figures  4-3,  4-5)   indicate  that  the 
major  error  of  the  model  is  the  reference  (hs)   fluid  whose 
composition  dependence  is  too  strong,  rather  than  too  weak, 
particularly  at  compositions  close  to  pure  component.  As 
a  result,  improvements  in  modelling  must  seek  a  funda- 
mentally different  concept  from  the  rigid  body  reference 
idea . 


We  now  examine  classical  expressions  for  activity 
coefficients  in  search  of  guidance  in  the  modelling  of 
direct  correlation  function  integrals.     We  analyze  a  few 
semiempirical  equations:     Van  Laar   (Prausnitz,  1969), 
Solubility  Parameter  (Prausnitz,   1969) ,  Wilson  (Prausnitz, 
1969),  NRTL   (Prausnitz,   1969),   and  UNIQUAC   (Abrams  and 
Prausnitz,  1975) .     We  pay  little  attention  to  the  physical 


6 . 4     Suggestion  from  Excess  Gibbs  Energy 
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meaning  of  the  parameters  involved;  our  main  interest  is 
in  the  compositional  form  of  the  equations . 

The  activity  coefficient  expression  from  the  Van  Laar 
model  is 


In  =   2  (6-12a) 

Taking  the  derivative  of  equation  (6-12a)  with  respect  to 
mole  fraction  of  component  1  and  introducing  this  result 
in  equations   (3-6)   to   (3-8)   we  get 

—2  2 

A'B' 
[X2+B'x^] 


V,V„         XX-  A'B' 

'-  ^  '  ^^^^ 

V-  x^  A'B' 


where  A'  =  2A  and  B '  =  2B 

In  solubility  parameter  theory,   the  activity  coeffi- 
cient has  the  following  expression 

-       2  2 
V,  x_  v„  (6,-6_) 
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where  v?  is  the  molar  volume  of  component  i  and  6j_  is  the 
solubility  parameter  which  is  calculated  using 


V 

6.   =  ^_Jh_ 
o 
^1 

with  AU^^  as  the  residual  internal  energy.  Rearranging 
equation   (6-16)   and  regrouping  terms,  we  can  recover 
equation  (6-12a) 


In        =  J  (6-12b) 

A  "^1 

with  A  =  v°   (<S;l"'^2^  ^  "  ^l/'^^ 

Because  equation  (6-12b)   is  a  particular  form  of 
equation  (6-12a) ,  the  direct  correlation  function  integrals 
for  the  solubility  parameter  theory  are  given  by  equations 
(6-13)   to   (6-15) .     A  similar  result  is  obtained  from  the 
three  suffix  Margules  equation.    All  of  these  are  derived 
from  the  Wohl  expansion  (Prausnitz,   1969) . 

In  the  Wilson  model,  the  activity  coefficient  expres- 
sion is 


In  y.  =  -  ln(xT+x,  G,-)   +  x-[ 


1  1     2  -12'       ^2^X3_+G,2X2  "  ^IV^ 


(6-17) 
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Following  the  same  procedure  as  for  the  Van  Laar  model,  we 
get  the  equations  for  C^^^,  C^2'  '^22' 

A  2         C.  X, 

(1-C.,)=—: ^+   ^  j+     Z   — =  J  (6-18) 

11         vkRT        (x^+x^B)^       i=l  (x^+X2D.)2 


V,V-         x,x-  A  2  C.  X, 

1^  VKKi  (X   +x   B)^         i=l    (X   +X  D.)^ 


V?  x^  A  2  C.  X- 

22         VKRT        (x^+x^B)^       i=l  (X^+X2D.)^ 


2 

where  A  =   (1-G^2)    '  ^  =  ^12'  ^1  ~  ^1  "  ^'  ^2  ~  ^^'^ 

°2  =  1/^1 

The  activity  coefficient  expression  from  the  NRTL 
model  is 

In  Yi  =         [  ^1     ^  +   1^     1^  2  1  (6-21) 

(x^+X2G2i)  (X2+X^G3_2) 

Taking  derivative  with  respect  to  mole  fraction  of  component 
1  and  using  equations   (3-6)   to  (3-8) 

2        A.  x^ 
1=1  (x^+B^x^) 


V,V-         2     A.  XX 

1=1 (x^+B^X2) 
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Where        =  t^^G^^^,         =  G^^,         =  t G^^^ ,  and  B2  =  G^^ 

In  the  UNIQUAC  model,   the  activity  coefficient 
equation  is 


m       =  m       +  (|)  qi  in  ^  +  *2  ^^1  " 


T  T 

^1  ^1    l^pi  +  o'^  e'+9'T       ^  (^-25) 

1     1     9i+92  21  ®2^^1^12 

The  direct  correlation  function  integrals  that  we  get 
from  equation  (6-25)  are 

4        X2A^  2  A. 

(1-c^^)  =       +  z   (X-+B.X-)         ^2     ~  ri 

1=1       1     1  2         j=l  (x^+BjX2) 


+  X2  J  (6-26) 

(x^+DX2) 


V,V~         4         x^A.  2  A. 

(1-C     )    =     ^  2  _     V  J-  J-         _     7     V  J 

"-12^        vkRT        .    ,    (x,+B.x_)  ^1    ,     ~  ~J 

1=1       1     1  2'        3=1  (x  +B.x_) 

X     3  z 


+  x^X2 


(x^+DX2) 


(6-27) 
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(I-C22)  = 


+  X 


2  C 


(6-28) 


where  the  parameters  of  equations   (6-26)   to   (6-28)   are  a 
complex  combination  of  the  parameters  of  equation   (6-25)  . 

The  Wilson,  NRTL,  and  UNIQUAC  are  local  composition 
models.    These  models  describe  excess  properties  only.  For 
representing  the  thermodynamic  properties  of  mixtures  at 
high  pressure,   it  is  convenient  to  use  an  equation  of  state 
in  all  the  ranges  of  densities.     With  the  equation  of  state 
approach  the  election  of  references  is  avoided,  and  mixtures 
containing  supercritical  components  can  be  described  without 
the  use  of  hypothetical  reference  state.     The  idea  of 
incorporating  the  local  composition  concept  into  an  equa- 
tion of  state  has  recently  been  the  subject  of  study 
(Whiting  and  Prausnitz,   1982;  Mollerup,  1981)  . 

An  equation  for  the  excess  Gibbs  energy  from  such 
equations  of  state  is 


(6-29) 

Whiting  and  Prausnitz  (1982)  developed  equation  (6-29)  as  an 
example  of  the  procedure  of  including  the  local  composition 


g     =   (-RT/a)      E     x.    In     Z     x.   exp   [ a (a . . -a . . ) /bRT] 
i=l     ^         j=l  ^ 
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vkRT 

(x^+X2B) ^ 

^1^2 

x^x^  A 

v<RT 

(x^+X2B) ^ 

^1  . 

x^  A 

concept  into  an  equation  of  state.     The  applicability  of 
equation  (6-29)   is  very  limited  because  it  was  assumed 
that  the  molecules  of  the  mixtures  were  of  approximately 
the  same  size. 

Using  the  second  derivative  of  equation  (6-29)  with 
respect  to  component  1  in  equations   (3-6)   to  (3-8) 

—2  2 

V,  X    A  2         C.  X- 

(1-^11^   =  ^      —772       .\         '     '     2  (6-30) 

  1=1  (x^+X2D^) 


<^-=12'   =  TjRf  -   ,     :  ^,2  -  X    ,  ,2  (6-31) 

  1=1  (x^+X2D^) 


(^-^2)   =  ^  ^   ,\     ,,2  +  ,  ,2  (6-32) 

(x^+X2B)         1=1  (x^+X2D^) 


with  A  =  (1-exp  [a(a2^-a^j_)/bRT] )  ^/a,  B  =  (exp 
[a(a2^-a^^)/bRT])/a,        =  -1/a,        =  B,  C2  =  C^,  and  D2  = 
(exp  [a(a^j^-aj_2)/bRT]  )/a 

It  is  evident  that  equations   (6-30)   to   (6-32)  are 
identical  to  equations   (6-18)   to   (6-20) .     Then,  the  form 
of  the  direct  correlation  function  integrals  for  this 
particular  local-composition  equation  of  state  is  the  same 
as  those  for  the  Wilson  model. 

We  do  not  pursue  further  in  the  search  of  forms  with 
different  and  more  complex  equations  of  state  because  we 
assume  they  may  produce  direct  correlation  function  inte- 
grals similar  to  those  from  Wilson,  NRTL,  and  UNIQUAC. 
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It  is  well  known  that  the  Wohl  form  and  the  Wilson 
form  represent  extremes  of  compositional  behavior  of  In  . 
The  former  shows  liquid-liquid  behavior  with  relatively 

oo 

small  values  of  In  y^,  whereas  the  Wilson  equation  cannot 
yield  immiscibility .     The  other  local  composition  models 
are  in  between. 

Comparison  of  the  appropriate  equations  in  (6-13)  to 
(6-32)   show  that  the  trend  for  C^^  is  from  a  variation  like 

A..  x2 

-iJ  :i 

( l+B .  . X  . ) ^ 

ij  : 

to  terms  of  the  form 

A..  x" 
(l+B. .   X.) ^ 

i:n  j' 

where  j  7^  i  and  n  =  1,2.     Clearly  these  all  include  binary 
parameters  for  C^^  (x^  ^  1) .    The  major  change  in  the 
representation  is  to  include  a  less  strong  (quadratic 
rather  than  cubic) ,  but  more  flexible   (more  terms)  varia- 
tion of  C . .  with  X . . 

Using  these  ideas,  we  can  suggest  an  alternative 
approach  in  modelling  C^^ .     As  it  is  very  convenient  to 
model  small  numbers  and  the  term  that  contains  the  partial 
molar  volvime  and  the  isotheinnal  compressibility  is  numeri- 
cally the  most  important,  we  suggest  to  model  the  difference 

(1-C..)   -   (1-C..)^^.     The  properties  present  in  (1-C..)^'^ 
•J  -'-J  li 
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are  pure  component  properties  that  make  the  computation  of 
the  ideal  direct  correlation  function  integrals  easy.  For 
modelling  this  difference  we  suggest  the  following  form 

AC^,  =  -    „    I        =     "      „       "j'  (6-33) 

RT  U^+B^^x^) 

The  equations  for  ^C^2  ^^'^  ^^22  similar  to  equation 

(6-33) . 

Finally,  we  must  point  out  two  serious  weaknesses  of 
this  suggestion.     First,  it  lacks  a  theoretical  basis. 
The  parameters  have  no  physical  meaning,  and  they  need  to 
be  evaluated  from  binary  experimental  data.     Second,  the 
extension  to  multicomponent  systems  is  highly  uncertain. 
On  the  other  hand,  the  main  advantage  is  that  the  calcu- 
lated properties  are  not  very  sensitive  to  the  particular 
form  of  the  direct  correlation  function  integrals  and  the 
results  can  be  accurate. 


CHAPTER  7 
CONCLUSIONS 


In  this  work  we  have  used  a  unique  approach  to  calcu- 
late liquid-phase  thermodynamic  properties.     The  fluctuation 
solution  theory,  as  we  call  it,  derived  from  rigorous 
statistical  mechanical  relations,  expresses  derivative 
properties  as  functions  of  direct  correlation  function 
integrals . 

The  thermodynamic  properties  of  interest  (activity 
coefficients  and  pressure)  are  obtained  by  integrating  these 
relations .     Other  approaches  get  the  same  properties  by  dif- 
ferentiation of  excess  Gibbs  energy  models  or  P-v-T-x  equa- 
tions of  state. 

The  most  important  assumption  made  in  fluctuation 
solution  theory  calculations  is  that  the  direct  correlation 
function  integrals  may  be  modeled  with  equations  of  simple 
form.     The  foundation  of  such  assumption  is  that  it  has 
been  found  that  for  dense  fluids  these  integrals  are  ap- 
parently insensitive  to  the  details  of  the  intermolecular 
interactions   (Gubbins  and  O'Connell,   1974)  . 

The  Mathias  Van  der  Waals  model   (1978)   has  been  ex- 
tended for  use  in  describing  the  solubility  of  gases  in 
binary  solvents  and  solubility  of  liquids  in  liquids .  We 
have  made  modifications  in  the  model  to  improve  the 
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modelling  of  liquid-liquid  equilibria.     In  the  process  we 
have  made  an  extensive  determination  of  direct  correlation 
function  integrals  from  experimental  data  to  find  out  how 
these  integrals  really  do  vary  with  composition.  Finally, 
we  have  made  some  suggestions  to  improve  the  model  for 
further  work. 

In  spite  of  the  relative  lack  of  success  in  describing 
the  liquid-liquid  equilibria,   the  approach  has  considerable 
potential.     The  assumption  of  simplicity  still  is  valid,  and 
it  is  clear  that  the  direct  correlation  function  integrals 
are  insensitive  to  the  details  of  the  intermolecular  forces. 
Whatever  complexity  appears  with  some  of  the  experimental 
direct  correlation  function  integrals  is  probably  due  to 
uncertainty  in  the  correlating  experimental  data. 

The  results  of  Chapter  4  that  extend  early  calculations 
of  Mathias  prove  that  the  model  is  excellent  in  describing 
gas-liquid  equilibria.     Unfortunately,  the  application  of 
the  model  in  liquid-liquid  calculations  shows  some  weak- 
nesses.    For  example,   the  perturbed  hard-sphere  model  does 
not  have  precisely  the  right  temperature  dependence.  Thus, 
closed-loop  systems  cannot  be  correlated  even  with  parame- 
ters linearly  dependent  on  temperature . 

Results  for  ternary  systems  indicate  that  the  composi- 
tion dependence  in  the  model  is  inadequate.     The  main  dif- 
ficulty seems  to  be  that  the  model  uses  only  one  binary 
parameter  whereas  more  are  necessary. 
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Further,  at  infinite  dilution  in  binary  systems,  data 
show  that  the  DCFI  strongly  depend  on  interactions  between 
the  components,  but  the  model  has  no  term  to  account  for 
these  interactions . 

A  suggestion  is  made  to  overcome  the  infinite  dilution 
problem.  It  may  also  apply  to  multicomponent  liquid-liquid 
equilibria. 

The  general  concept  is  to  incorporate  an  additional 
term  in  the  expansion  of  the  perturbed  hard-sphere  model. 
The  term  is  quadratic  in  composition  and  resembles  the 
third  term  of  the  virial  equation  of  state.     With  this 
modification  the  interaction  between  components  at  infinite 
dilution  is  incorporated  into  the  model,  and  the  necessary 
binary  parameters  are  available  to  describe  liquid-liquid 
equilibria . 

Another  suggestion  for  modelling  direct  correlation 
function  integrals  is  to  use  the  general  forms  obtained 
from  excess  Gibbs  energies.     In  such  an  expression,  the 
parameters  have  no  physical  meaning  and  they  need  to  be 
evaluated  from  experimental  data.    As  a  result,  this  is  less 
satisfying  than  the  above  idea. 

It  is  hoped  that  the  present  investigation  will  have 
provided  a  basis  for  analysis  of  data  and  establishment  of 
models  of  liquid  solutions  that  will  lead  to  more  accurate 
correlations  and  predictions  of  these  most  complex  systems. 


APPENDIX  A 
EQUATIONS  IN  THE  BASIC  MODEL 


Introducing  (4-6)  in  the  equations  of  Chapter  2,  we 
obtain  the  thermodynamic  properties  for  a  multicomponent 
fluid  phase. 

The  pressure  equation  is 


=  P^^T,pf)    -  P^^T,p^) 


RT 


where  the  hard  sphere  term  is  the  Carnahan-S tar ling  (1969) 
equation  for  a  hard  sphere  mixture 

V  (T,p)   =  I  [       °y  +       ^  ^  3  +   1-  ]  (A-2) 

^     ^     ^3^        (1-^3)  (1-^3) 

where 


=  f  Pia^  k  =  0,1,2,3  (A-3) 


and  where  f^^,  g^^  are  defined  in   (4-9a) ,    (4-9b) ,    (4-9c)  . 

The  activity  coefficients  have  the  following  expres- 
sion 
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In 


f        hs   ,„    fs         hs  r> 


r  RT  RT 


+  2        [p^-p^]v*j  (g^j-fij)  (A-4) 


where  the  hard  sphere  chemical  potential  is  given  by  the 
Carnahan-Starling  expression 


hs  ^  2  ^ 

^i  ^2^i  ?2^i 

(T,£)   =  In  p  -  In   (1  -  5  )  [1  -  (  )    [3-2  {-^)]] 


+   [3(?2^i  +  5l^i)  -  (52^i^^^2  -53)/C3]/(l-?3) 


*  * 

The  cross  characteristic  parameters,  V^^  and  T^^  are  defined 
by 


*  1/3        ^  1/3  ^ 

*  (V. .  +  V. .  ) 

*  *       *  1/2 

=   (1  -k^j)    (T^  -T.)  (A-7) 


and  the  function  f  and  g  are  given  by  the  following  expres- 
sion 
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.       0.65386227  „     *  ^  ^  , 

^  =   ,  0  .16067976  ^/^i  ^  ^'^^  (A-8) 

(T/T^) 


=  0.807662393  exp   [-0.22010926    (T/T*) ] 


T/T*  <  0.73  (A-9) 


g  =  0  .3625065  -  ^ '^^^J^^^  t*  .  -  1-75^3882  ^*  ^ 

T  ID  ip-  ij 


,1,470750  0  .0041793 


The  relation  between  and  that  is  needed  to  obtain  the 
hard  sphere  diameter  from  f  is 

„  3 

2n  ^i 

-3-  -S         =  fi  (A-ll) 

V.  .  ^ 
11 


APPENDIX  B 
HARD  SPHERE  EQUATIONS  FOR  MIXTURES 


The  Carnahan-Starling   (1969)   equation  was  used  by 
Mathias   (1978)   in  modelling  direct  correlation  function 
integrals.     However,  introducing  a  parameter,   a,   it  is 
possible  to  generate  other  equations  for  hard  sphere 
fluids .     The  pressure  equation  is 

hs  _  6kT  ^o  ^  ^1^2  ^2 

^3^2 

+  a   \  )  (B-1) 

where 


1  „  k 


5k  =  6  n     E     p^a/^  k  =  0,1,2,3  (B-2) 

i=l 


k  is  the  Boltzmann  constant,   p^  the  number  density,  and 

the  hard  sphere  diameter  of  component  i.     In  the  follow- 
ing equations,  the  superscript  hs  will  be  dropped,  but  all 
the  thermodynamic  properties  will  be  for  hard  spheres. 

The  pressure  and  chemical  potential  are  related 
through  the  Helmholtz  free  energy  by  thermodynamic  relations 


190 


191 


P  =  - 


9A 
3V 


(B-3) 


T,N 


9A_ 


{B-4) 


T,N. 


We  obtain  the  Helmholtz  free  energy  by  integration  of 
equation  (B-3) ,  and  then  we  use  equation  (B-4)   to  get  an 
expression  for  the  chemical  potential  of  a  hard  sphere 
mixture . 


3  3 
yj_      In        Aj_      n  P 

^  ^      ^int,i  " 


-In  (1  -  ^3) 


3  ?2 


3  Cia. 


9  ^2  a. 


^  (i-?3)    -nrriF^  — 
2 


^2 

?2 


-    (-rr^)  [in(i-?3)  j^-^ 


2(l-5o) 


C  a.  2 

a  (^-^)      [2  In   (1  -  C-,)   +  C 


(2  -  53) 
3   (1  -  ?,) 


]  (B-5) 


where  is  the  thermal  De  Broglie  wavelength,  and  q^^^^^  ^ 
is  the  molecular  internal  partition  function. 

The  relation  between  the  direct  correlation  function 
integrals  is  given  by  equation  (2-23) 


6.  .       C.  . 


(2-23) 


T,  p. 
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Using  equations  (2-23)   and  (B-5)  we  get  the  i-j  direct 
correlation  function  integral  for  a  multicomponent  hard 
sphere  mixture . 

-r    ^=  -  n  (      1        l/6na/  ,   3n  ^i 
^ij     -  P  ^  TTT^         3       T  TTip- 

.    311  ^j^^2^i  ^  311  '^j^i^     ^  3n  ^j^^l^i^ 
'    (I-C3)'        '  6  (W3)2 

^  2^  2  3-2  2 

T,  a.  a.  „  a.  Co  C7- 

+  95,  I    J     ^  ^  +  9  f    J  ^ 
^  ^   (I-C3)  ^  (1-^3) 

or  2  3- 

^3  ^3^ 
[ln(l-C  )   +  j^-L.y  1-^  ] 

^3^  2(1-^3)^ 
3 

^2*^1  ^     n         ^2^i       IT  ^ 


1 


h         ^       C32  6  ^3 

(2-C-,) 

(2  m  (1-C3)  +  C3  -n:r^  ) 

^2*^1   3  ^       <2-C-)  a.^ 

+  a  {-^)      [  n/6  a/c^   ^  -  2n/6  -rf- 

^3  =>     ^  (1-C,) 
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+  a  34  (i.r  )   +    °     J     ^    H-  3n/6 

^  6(1-^3)^  (1-^3) 

2  3_  3-   _  3 
a.  a.                    a.  ^^CjC^i 

+  3II/6  -3  ± — ^  +  6n/6   ^  "^^^ — 

(I-C3)  (1-^3)^ 

-  2     2     3  3^  3  3 

t,^  o  .  a .  a  .  E,y  a  . 

+  9n/6  -= — 2 — i_  +  9n/6  — = — J — 

(1-C3)  (1-C3) 

3     _  3  3     3_   _  2 

a .  a  .  Co  a .  a  .  C-^Co 

+  a[  n/6  — — +  3n/6  — — J  -^^"^ 
(l-?3)  (I-C3) 

3  3-  -  3 
a.  a.  53^2 

+  3n/6  -J — i — f-^  ]}  (B-6) 
(I-C3) 


For  a  pure  component  system  equation  (B-1)   reduces  to 


p  =  kTp  (1     n  +  n^+an^)  (b.7) 
(i-n) 


where  x],   the  packing  fraction,   is  ^3-     Equation   (B-7)  is 
equation  (4-14)   in  Section  4.3. 

In  equation  (B-7)  when  a  is  equal  to  -1,   the  Carnahan- 
Starling  equation  is  recovered,  when  a  =  0  the  Percus- 
Yevick  compressibility  equation  is  recovered,  and  when 
a  =  -3  the  Percus-Yevick  virial  equation  is  recovered. 

The  direct  correlation  function  integral  for  a  pure 
component  is 


c  =  -8n  +  2n^  -  4(i+a)       +  (i+a)  (g.g^ 
(1-n)^ 
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In  the  particular  case  of  Carnahan-S tar ling,  equation  (B-8) 
becomes 


c  = 


2n(n-4) 


(i-n)"^ 


(B-9) 


APPENDIX  C 

GIBBS-DUHEM  ANALYSIS  OF  THE  MULTIFLUID  APPROACH 

The  starting  point  in  the  multifluid  approach  is  equa- 
tion (5-1) 

n 

In  Y     =     E     X    In  ^y.  (C-1) 
1       j=l     3  1 

where  -'y^  is  the  activity  coefficient  of  component  i  using 
reference  j ,  taking  derivatives  on  both  sides  of  equation 
(C-1)  at  constant  T  and  p 

n  .  n  . 

d  In  y.  =     I    X.  d  In  ^Y-  +     ^     d  x.   In  ^y.  (C-2) 
^      j=l     3  ^      j=l        3  ^ 

If  the  Gibbs-Duhem  equation  is  applied  to  the  activity 
coefficients  defined  by  equation  (C-1) ,  we  can  write  using 
equation  (C-2) 


n 

I  X.  d  In  Y-  =  0 
i=l  ^ 


n  n 

Z  X.  Z  X.  d In  ^Y4 
i=l        j=l.  ^ 


n  n  . 

+     E    X.     E     dx.  In-'y- 
i=l     ^  j=l        3  ^ 


(C-3) 
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As 


n 

Z  X 
i=l 


n 

Z 


.  d  In  ^■ 


X.  a  m 


n  n 

Z    Xj  Z 


j=l 


i=l 


X .  d  In    Y  • 
1  'i 


=  0 


(C-4) 


because 


n 

Z  X.  d  In  \.  =  0  (C-5) 
i=l    ^  ^ 


by  model  construction,  then  for  equation  (C-3)   to  be  true 


n  n 

(C-6) 


Z  X.  Z  d  X.  In  -"y.  = 
i=l     ^  j=l  3 


as  the  mole  fractions  are  not  independent,  equation  (C-6) 
can  be  rearranged 


n  .  , 

Z    X.     Z    d  X.    (In  -  In     y.  ) 

i=l    ^  ijik        ^  ^  ^ 


n  . 

Z  d  X.  (  Z  X.  (In  ^Y.  -  In  y,) 
jfik        ^     i=l    ^  ^  ^ 


=     0  (C-7) 


In  equation  (C-7)  the  mole  fraction  differentials  (d  x ^ ; 
jT^k)  are  independent,  then 
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n  .  , 

Z  X.  (In  ^Y,-  -  In  yJ  =0  for  all  j^^k  (C-8) 
i=l    ^  1  1 


Forming  a  linear  combination  of  (C-8)  equations 


n 

X  . 


Z    X.     Z    X.    (In  -"y.   -  In    Y-)   =  0 
^  i=l     ^  1  1 


or 


n  n  ^  n  , 

In    Yi  (C-9) 


Z    X.     Z    X.   1     -"y-   =  Z 
j=l     3  i=l    ^    ^  i=l 


or 


n 

-'y^     for  any  k  (C-10) 


In    Yi  =       2    X.  In  -"y. 
1        i=l     3  1 


Equation  (C-10)   says  that  the  only  possible  choice  of  ref- 
erence that  will  satisfy  the  Gibbs-Duhem  equation  is  a  pure 
component  reference. 

We  can  extract  the  same  conclusion  by  starting  with 
the  definition  of  excess  free  Gibbs  energy  in  the  multi- 
fluid  approach 


gE 

RT  =  ^i   1^  ^i 


n  n 

=     Z    X.     Z    x     In  ^Y.  (C-11) 
i=l     1  -1=1     3  1 
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By  definition  of  the  activity  coefficient  (which  assximes 
the  Gibbs-Duhem  equation) 


In  Yi 


(  ^  ) 


3N,     '  RT 


P,T,N 


(C-12) 


where        are  the  number  of  moles  of  component  i  and  N  the 
total  number  of  moles,  using  equation  (C-11)   in  equation 
(C-12) 


As 


3     (  N_2l  ) 


8Nj^   '  RT 


P,T,N 


P,T,N 


n 


n 


n 


?  i=l     ^  j=l     3  ^1      N   .^^  : 


n 


-,  n  h  9  In  -'Y• 
+  ^     E     N.      Z     N.   rr^  ^ 


(C-13) 


P,T,N 


n 


n 


-EN.  EN. 
N  i=l     ^  j=l  3 


3  In  ^Y. 
3N 


1  ^ 

=  77     E     N.  E  N. 

N    .    ,      1  .  1 

3=1  1 


k 

3  In  ^Y. 

3N. 
1 


=  0 


(C-14) 
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Equation  (C-13)   reduces  to 

n  .  n  k  "  i 

In  Yj^  =     2     X  .   In  -"y    +     Z  [  In         -  (  E     x  .  In  -"y.  )  ] 

j=l  i=l  j=l 

(C-15) 

but 


n 

S  X .  In  -'yi  =  In  Yi 
j=l     ^  ^  ^ 


Then,  we  can  write  equation  (C-15)  as 


^  k  "  -i 

In  Yv  =  In  Yv  +     2     x.    [In    y.  -  (  ^     x.  In  ^y. )  ]  {C-16) 
^  ^      i=l  ^       j=i     3  1 


Equation  (C-16)  can  be  true  if 


n  .  n  . 

E    X    [In  -   (  Z    X.   In  ^y . )  ]   =  0  (C-17) 

i=l    ^  1        j=i     J  1 


or 


k  "  i 

In    Y^  =     2     X.  In  -"y.         for  any  k  (C-18) 
1      j=l     3  1 


Equation  (C-18)   is  the  same  as  equation  (C-10) 
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